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ABSTRACT 


We wish to solve fluid flow problems in only a portion of a large or infini te 
domain. By restricting our area of interest, we effectively create a boundary where 
none exists physically, dividing our computational domain from the rest of the physical 
domain. The challenge we must overcome, then, is defining this boundary in such a 
way that it behaves computationally as if there were no physical boundary. Such a 
boundary definition is often called a non-reflecting boundary, as its primary function 
is to permit wave phenomena to pass through the open boundary without reflection. 
In this dissertation we develop several non-reflecting boundary conditions for the 
linearized Euler equations of inviscid gas dynamics. These boundary conditions are 
derived from the Higdon, Givoli-Neta, and Hagstrom-Warburton boundary schemes 
for scalar equations, and they are adapted here for a system of first-order partial 
differential equations. Using finite difference methods, we apply the various boundary 
schemes to the gas dynamic equations in two dimensions, in an open domain with 
and without the influence of gravity or Coriolis forces. These new methods provide 
significantly greater accuracy than the classic Sommerfeld radiation condition with 
only a modest increase to the computation time. 
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I. INTRODUCTION 

Many problems in compntational flnid dynamics occnr within a limited portion 
of a very large or infinite domain. Difficnlties immediately arise when one attempts 
to dehne the bonndary condition for snch a system. Snch bonndary conditions are 
necessary for the problem to be well-posed, bnt the physical system nnder consider¬ 
ation has no bonndary to model. One needs to dehne an artihcial bonndary whose 
behavior models the open edge of the physical system. Snch a bonndary dehnition is 
often called a non-reflecting bonndary condition (NRBC), as its primary fnnction is 
to permit wave phenomena to pass throngh the open bonndary withont rehection. 

NRBC development has been an ongoing research area since the 1960s. As with 
any compntational endeavor, there are trade-offs involved. An ideal NRBC wonld be 
fast, accnrate, stable, and easy to implement: fast, meaning that the compntation 
of bonndary valnes is small or negligible relative to the interior domain; accurate, 
meaning that there is little to no spnrions rehection indnced by the bonndary con¬ 
dition; stable, meaning that the bonndary compntations do not canse the solntion 
to degrade catastrophically over time; and easy to implement, meaning that the nser 
can incorporate the NRBC compntations into an operational model with minimal 
modihcation to existing code. Realistically, one mnst settle for two or three of these 
attribntes, at best. Conseqnently, the search for better NRBC methods continnes. In 
addition, researchers continne to apply existing NRBC methods to new domains and 
wave propagation eqnations. 

In this dissertation we develop several NRBCs for the linearized Enler eqna¬ 
tions of inviscid gas dynamics. These bonndary conditions are derived from the 
Higdon, Givoli-Neta, and Hagstrom-Warbnrton bonndary schemes for scalar eqna¬ 
tions, adapted here for a system of hrst-order partial differential eqnations. Using 
hnite diherence methods, we apply the varions bonndary schemes to the gas dynamic 
eqnations in two dimensions, in an open domain with and withont the inflnence of 
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gravity or Coriolis forces. These new methods provide signihcantly greater accuracy 
than the classic Sommerfeld radiation condition with only a modest increase to the 
computational cost. 

Our motivation for developing these NRBCs is to support the work of Giraldo 
and Restelli in their efforts to develop the next generation of mesoscale atmospheric 
modeling tools [32]. In their current form, the models rely on large absorbing layers 
surrounding the computational domain. In order to be effective, these “sponge lay¬ 
ers” can be as thick as the original domain [31], resulting in a total domain which 
is nearly quadrupled in size. If we can replace these large sponge layers with accu¬ 
rate NRBCs, then the modeling tools will requires less memory and computational 
overhead, signihcantly increasing their efficiency. 

The rest of the dissertation is outlined as follows. We begin in Chapter II by 
deriving the equations under consideration, the linearized Euler equations of inviscid 
gas dynamics. In Chapter III we provide an overview of existing NRBCs and provide 
specihc details about the Higdon, Givoli-Neta, and Hagstrom-Warburton schemes for 
scalar equations. We then extend these boundary conditions to the hrst-order 2-D 
linearized Euler system in Chapters IV (Higdon), V (Givoli-Neta), and VI (Hagstrom- 
Warburton). In all three cases, we consider the NRBC implementation in a semi- 
inhnite or inhnite channel and in an open domain, under basic conditions as well 
as under the inhuence of Coriolis forces or gravity. Numerical examples using hnite 
differences are provided throughout. We discuss the issue of long-time stability in 
Chapter VH. We offer a qualitative comparison of the three NRBC techniques in 
Chapter VHI. We close in Chapter IX with a summary of our results and a list of 
areas for further research. 
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II. 


MODELING INVISCID FLUID FLOW 


In this chapter we explore the equations governing the motion of a body of 
fluid. These principles describe the flow of water in a channel or in the ocean, air 
movement over mountains, airflow and drag in aircraft design, and even the heat 
generated by a spacecraft re-entering the atmosphere. Although many physical phe¬ 
nomena depend on the viscosity of the fluid, certain large-scale flows of low-viscosity 
fluids (e.g., air) can be reasonably approximated by assuming the viscosity is negli¬ 
gible. 

By neglecting viscous forces, our fluid flow equations can be derived based 
simply on physical conservation laws governing mass, momentum, and energy. The 
following section considers each conservation law in turn and derives the relevant 
governing equations therefrom. 

We derive the Euler equations based hrst on internal factors. Then we consider 
the inhomogeneous factors which affect these equations in the context of atmospheric 
modeling. 


A. DERIVATION OF EQUATIONS 

1. Conservation of Mass 


Mass is conserved in a closed system. If dm denotes an inhnitesimal portion 
of the mass, then / dm denotes the total mass of the body, and conservation of mass 
requires 


d 

dt 


dm = 0 


(II.l) 


Since the mass is continuous, we can bring the derivative inside the integral; thus 



( 11 . 2 ) 
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Furthermore, since this statement must be true for any piece of the body’s mass, the 
integrand must be identically zero, i.e.. 


or 


— (dm) = 0 
dV ^ 


— {pdxdydz) = 0 

CLL 


Applying the product rule gives 


/ 7 7 7 \dp / 7 7 \ / 7 7 \ AV / 7 7 \ A^ 

(dxdydz)— + (pdydz)d— + (pdxdz)d— + (pdxdy)d— = 0 

Dehne the velocity vector u = (u, n, where u = dx/dt, v = dy/dt, and w 
and divide all four terms by dxdydz to get 

dp du dv dw 

Applying the chain rule to each derivative gives 

7^+U^+V^+wlf+pip'll + + ' 

at ox oy oz • \ox ax oy dx oz dx J 

-Ln (\ ^ (I dtv^ i dw dz \ 

\dx dy dy dy ' dz dy) ' r \ dx dz ' dy dz ' dz dz) ) 

However, since x, y, and ^ are independent, this equation reduces to 
dtp + udxp + vdyp + wdzp + p {dxU + dyV + d^w) = 0, 


y =0 


i.e., 


or, in vector notation. 


&tP + &x{pu) + dyipv) ^ d^ipw) = 0 


ftp + V ■ (pm) = 0 , 


(11.3) 

(11.4) 

(11.5) 
dz/dt, 

(11.6) 

(11.7) 

(11.8) 

(11,9) 

( 11 , 10 ) 


where from here on, we use the following shorthand for partial derivatives: 

f, _8_ a 

* dt ' dxdy 

This equation assumes no external sources or sinks adding or removing mass from 
the system. We note in passing that we have shown, from first principles, a special 
case of Reynold’s transport theorem [3] applied to the density of a fluid. 
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2. Conservation of Momentum 

We now consider the momentnm of the flnid body. To this end, we hrst 
determine what forces act npon the flnid. For the pnrpose of this dissertation, the 
flnid body is assnmed to be a portion of the Earth’s atmosphere. 

a. Forces Acting upon a Fluid Body 

The internal forces exerted by the body on itself consist of pressnre 
forces and viscons forces (which are neglected since we are assnming an invsicid 
flnid). The external inhomogeneons forces we consider are the effects of gravity and 
the Earth’s rotation. 

(1) Force Dne to Internal Pressnre. The pressnre force 
resnlts from pressnre differences within the body (see Fig. 1) and acts to retard the 
motion of the flnid. 



Fignre 1. Pressnre Differences in a Small Volnme [From [113], Fig. 3, p. 15] 

In a small piece of flnid volnme dV = dxdydz, the change in 
pressnre in the x direction is (to hrst order accnracy) 

dxP dx 
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The force exerted by this pressure difference is the product of the pressure difference 
and the surface area on which the pressure is exerted. Hence, the pressure force in 
the X direction is given by 


dFpressurex dx^ dydz ^xP dV 

Similarly, the forces exerted by pressure differences in the y and ^ directions are given 
by 

dFpressurey = “ {dyP dy) dxdz = -dyP dV 

dFpx^ggi^xcz i^^zP dz^ dxdy ^zP dV ^11.12^ 


Thus, the total force due to internal pressure differences is the sum of the three 
components. 


dFpressure = “ {dxP^ + dyPJ + dzPk) dV = -{Vp)dV, 


and the overall force on the entire body is 


pressure 


J dFpressure J {^p^dV 


(11.13) 


(11.14) 


(2) Force Due to Gravity. According to Vallado [112], the 
gravitational acceleration comes from Newton’s Law of Gravitation: 


a = 


/i 


r 

If I 


(11.16) 


where p is the Earth’s gravitational parameter (398,600.4418 km^/s^), and ||f|| is the 
radius from the center of the Earth (6378.137 km at sea level). At sea level, this gives 
a gravitational acceleration of approximately 9.798 m/s^. Gomparing the magnitudes 
of the acceleration at sea level to that at an altitude of 20 km, we find 


Aa = /i 


1 


1 


(r + 20km) 2 


m 

-0.06116 — 


(11.16) 


Given such a small relative difference, we can treat the gravitational acceleration as 
constant, so the gravitational force is simply 


gravity 


-gdmk 


(11.17) 
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where 


g ^ 9.798— (11.18) 

s 

(3) Force Due to Earth’s Rotation. In Fig. 2, D repre¬ 
sents the angular velocity of the sphere, which for Earth has a value of 27r radians per 
sidereal day (23 hours, 56 minutes, 4.090524 seconds [112]) or 7.292116 x 10“® rad/s. 
We must consider the effects caused by the rotating reference frame. Here we follow 



Figure 2. Rotating Sphere [From [113], Fig. 2, p. 11] 

the derivation laid out by Holton [70] and Pedlosky [94]. Let r denote the position 
vector of an element in our rotating reference frame (see Fig. 3, where 'ijj is the angle 
between Fand the angular rotation vector D). In a small time interval At, the vector 
rotates by an angle A6 = |D|At. If the time interval is small enough, the change in 
the vector is given by 

Af= n\f\A6 sin'ijj, 

where n is the unit vector in the direction of the cross product fl x f, i.e., 

^ D X r 

n = — - 

\nxr\ 


(11.19) 

( 11 . 20 ) 
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Figure 3. Vector in Rotating Reference Frame [After [94], Fig. 1.5.2, p. 15] 


Thus, 


or in the limit as At —)■ 0, 


Ar 


Q X f 
X r] 


rl-—sin-iw, 
' 'At 


df 

dt 


Q X f 
X r] 


|r1 


dO 


( 11 . 21 ) 


( 11 . 22 ) 


and since ^ = |fl|, we get 


dr ^ ^ 

— = Vt X r 
dt 


(11.23) 


If we have f = xi + yj + zk in the east-north-up coordinate frame shown in Fig. 2, 
then we convert from rotating to inertial reference frames by 


^ dr ^ ^ 

ui = u-\- — = u + \lxr^ 
dt 


(11.24) 


where u is the change in position within the rotating reference frame, and uj is the 
same motion relative to the inertial reference frame. Differentiation of the above. 
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substituting uj for r and noting that is constant, yields 


dui\ / dui\ - ^ 

(j ‘7 / —^ ^ ^ 

= — + 2nxu + nx{nxf^ 

CLL 


(11.25) 


Based on Fig. 2, our angular momentum vector is 


12 = 1121 cos (pj + 1121 sin (f)k 


(11.26) 


Using this dehnition and the dehnition of the cross product twice, we have 

/ a; \ 

n X (12 X r) = |12|^ ^ sin 0 cos 0 — 1/ sin^ 0 , (11.27) 

\ —^ cos^ (f) + ysmcf) cos (f) / 

but, since |12p = 5.317496 X 10 we can ignore the last term of (11.25) to get 


/ dui\ du 

-r- = — + 212 X n 

\ dt J j dt 

Applying this formula to the momentum vector gives 


(11.28) 


= ^{udm) + 212 X (udm) (11.29) 

CLL 

Hence, when we consider the change in momentum, we must add to it this rotational 
effect. Now this rotational effect can be simplihed. The atmosphere is thin compared 
to the radius of the Earth. Furthermore, atmospheric flows tend to be parallel to the 
ground (i.e., the velocity in the ^ direction is small). Taking the cross product of this 
vector and our velocity vector gives 


2(u,drn) 


212 X -u = 2|12|((wcos0 — usin0)£ + usiiKpj — u cos (j)k) 


(11.30) 


Again, assuming a thin atmosphere and thus neglecting terms in the ^ direction, this 
simplihes to 

212 X M 2|12| sm(p{—vi + uj) = 2|12| sm(j{k x u) (11.31) 
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If we define / = 2|fl| sin0, we can express this force as 


d 




dt 


{udm)\ + f{lv X (udm)) 




(11.32) 


R 


Since we are observing the flow from the rotating reference frame, this additional force 
becomes an inhomogeneons term on the right-hand-side of the momentnm eqnation: 


Fcorioiis = J f{ux k)dm 

b. Summary of Forces 

Thns, the eqnation for conservation of momentnm is 


(11.33) 


A 

dt 


J udm = — J WpdV -|- j f{u x k)dm — J 


gkdm 


Here, as before, the integrand mnst be identically zero. 


^iudm) -t- ^^dm = fiu x k)dm — gkdm, 
dt p 


(11.34) 


(11.35) 


or, by component. 


^(udm) + = fvdm 

dt p 

d din 

— {vdm) ^ - dyp = —fudm (11.36) 

d dm ^ 

— iwdm) - OzP = —gdm 

dt p 

Using a derivation similar to that in Section 1, these eqnations can be rewritten in 
vector notation: 


dt{pu) + V • {pu <^u) + Vp = f{pu X k) — gpk, 


(11.37) 


where 0 denotes the tensor prodnct. However, we can simplify the eqnations to 


dtu + udxU + vdyU + wdzU -\ — d^p = fv 

P 

dtv + udxV -|- vdyV + wdzV H —dyP = —fu 

dtw + udxW + vdyW + wdzW H —dzP = —g, 

P 


(11.38) 


or, in vector notation. 


dH 1 /V /V 

— -t- (m • V)m -f -Vp = fiu X k) — gk 

dt p 


(11.39) 
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3. Conservation of Energy 

The intrinsic energy in each piece of mass dm consists of its kinetic energy. For 
atmospheric modeling, there are also inhomogeneons components for thermal energy 
and potential energy. These qnantities are dehned by 

dHE = CyTdm 

IliTIP 

dKE = (11.40) 

dPE = gzdm 

In these eqnations, Cy is the specihc heat at constant volnme, and T is the temperatnre. 
Hence the total energy in the body is 

J {^yT + dm 

To change the amonnt of energy, apply force to the body over a distance or add/remove 
heat. Adding/removing heat is a complicated process involving solar radiation, ther¬ 
mal radiation from the Earth, heat dissipation into space, reflectivity of the Earth’s 
snrface, and other factors. Dne to the complexity, factors which contribnte to a change 
in temperatnre will not be considered here. Rather, we will simply consider the time 
derivative of onr heat energy term on its own withont determining the precise factors 
which determine it. Later in this section, we will hide the temperatnre component 
entirely, combining thermal and kinetic energy into a “total energy” term e. 

For force application, the momentnm eqnations assnme a balance between 
the external forces and the internal pressnre. However, if the two qnantities are not 
matched, then the volnme will expand/contract so that they are matched. If there 
is an imbalance between the forces on opposing sides of the body (hence a pressnre 
difference), then the body will accelerate in the direction of the net force. (Imagine 
pnlling a spring across a table. The spring will stretch according to the pnlling force, 
and the stretched spring will move in the pnlling direction.) 

We can conceive of the total force as having these two components, compression 
and acceleration. For compression at constant temperatnre, as the volnme decreases. 
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the internal pressnre increases, and vice versa. Hence, 




(11.41) 


However, if the temperatnre is not constant, then the pressnre increases proportionally 
to the thermal energy at the constant volnme, adding the following component to onr 


eqnations: 


^(c^Tdm) = ^dV 
dV ^ dt 


(11.42) 


The acceleration component is determined by Newton’s Second Law. The acceleration 
occnrs as a resnlt of an imbalance between the force acting on one side of the body 
and the force acting on the opposite side, which resnlts in the pressnre differences 
dehned in the conservation of momentnm section. Using F = ma and onr dehnition 


of dFpressure WC have 


-{'Vp)dV Vp 

dm p 


(11.43) 


This force changes the kinetic energy of the body. The time derivative of kinetic 


energy is 


— {KE) = dm iu ■ = dm{u ■ a) 


Using onr formnla for the acceleration, we have 


rj —nTTl 

— (iLE) =- {u ■Vp) = -{u- Vp)dV 

CLv p 


(11.44) 


(11.45) 


In the context of atmospheric modeling, this force also resnlts in an increase in the 


body’s potential energy as it changes position vertically, so that we have 

^^{KE) + j^{PE) = -{u-Vp)dV 


(11.46) 


Therefore, changing the energy in the body results in increased kinetic energy, body 
compression, and pressure increase due to temperature 


_ / ^ + sd dmU / -yv) = - /(«. Vp)dr + / liV (11.47) 
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Combining the integrals and bringing the derivative inside gives 


/1 + 

where we use the product rule to combine 
integral must be true for any body, so the 

d((c„r+l^+9.)dm) 


{u ■Vp)dV + p-^{dV) = 0, (11.48) 

^{pdV) — ^dV into p-^{dV). Again, this 
integrand must be identically zero 

+ {u-Vp)dV+ p^{dV) = Q (11.49) 

CLL 


II “'ll 2 

Let e = CvT + denote the total internal energy of the system (i.e., not including 
potential energy). Then we have 


A 

dt 


((e + gz)pdV) + {ud^p + vdyp + wdzp)dV + p^{dV) 

CLL 


0 


(11.50) 


It is easy to show that this equation is equivalent to 


dt{pe) + dx{{pe + p)u) + dy{{pe + p)v) + dz{{pe + p)w) = -gpw, (11.51) 
or, in vector notation, 

dt{pe) + V ■ ((pe + p)u) = -(gk) ■ (pu) (11.52) 

This equation can be placed in a simpler form, if we assume our fluid is an ideal 
gas. Using the ideal gas law p = pRT [24, 58], where R = Cp — c^, simple algebraic 
manipulation using the previous equations can reduce this energy equation to 


dtp + udxP + vdyp + wdzP + 7P {dxU + dyV + dzw) = 0 , (11.53) 

where 7 = —. In vector notation, this equation can be written as 

Cy 

dtp + {u ■ V)p + 7p(V • -u) = 0 (11.54) 

While this equation is simple, it lacks an explicit energy term and thus fails to convey 
the energy conservation principle from which it was derived. 
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B. SUMMARY OF NON-LINEAR EULER EQUATIONS 

Using only the physical conservation principles for mass, momentnm, and en¬ 
ergy, we can derive Enler’s eqnations for inviscid flnid motion. We have a system of 
hve eqnations for six variables: 

dtp + da;{pu) + dy{pv) d^{pw) = 0 

dt{pu) -h d^{pu^) + dy{puv) + d^{puw) + d^p = fpv 
dt{pv) + d^ipuv) + dy{pv‘^) + d^{pvw) + dyp = -fpu ( 11 . 55 ) 
dt{pw) + dx{puw) + dy{pvw) + dz{pw^) d^p = -gp 
dt{pe) + dxiipe + p)u) + dy{{pe + p)v) + dz{{pe + p)w) = -gpw 

In vector notation, the eqnations can be written as 

dtp + V ■ {pu) = 0 

dt{pu) + V ■ {pu ®u) + Vp = f{pu X k) — gpk (11.56) 

dt{pe) -f V ■ ((pe + p)u) = -{gk) ■ {pu) 

A state eqnation of some kind is reqnired to close the system. For an ideal gas, we 
can nse p = pRT to simplify the eqnations to 

dtp + d^{pu) + dy{pv) + d^{pw) = 0 

dtu -f udxU + vdyU + wdzU -\—dxP = fv 

P 

dtv + udxV -|- vdyV + wdzV H— dyP = —fu (11.57) 

dtw + udxW + vdyW + wdzW H —dzP = —g 

P 

dtp + udxP + vdyp -f wdzP + 7P {dxU + dyV + dzw) = 0 

(see Appendix A for details). For all three formnlations, the terms on the right-hand- 
side denote forces specihc to atmospheric modeling. 

In the next section, we will derive the linearized form of (11.57), which will 
form the basis for developing the hnite-difference implementation of the NRBCs in 
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Chapters IV-VI. We will begin with a simplified prototype implementation of the 
antomated Higdon NRBCs with no advection or forcing terms, and we will bnild 
the implementation to inclnde Coriolis, gravity, and non-zero mean flow. We will 
also develop anxiliary variable forms which eliminate the high-order derivative terms, 
nsing both the Givoli-Neta and the Hagstrom-Warbnrton anxiliary variable NRBC 
formnlations. 

C. LINEARIZED EULER EQUATIONS 

Having derived the non-linear eqnations, we now seek to simplify them into a 
linear form. We assnme that each state variable consists of a time-independent mean 
valne and a small pertnrbation from that mean. Thns, 

ip = (p + ip*, (11.58) 

where the overbar denotes the reference valne, and the asterisk denotes the 0(S) 
pertnrbation variable. Before we can derive this linearized form, we mnst first define 
onr reference variables. We will perform the linearization on (11.57). 

1. Defining the Reference Variables 

The reference valnes are time-independent by definition, bnt they may not 
necessarily be constant in space. It is reasonable to believe that the reference valnes 
for the velocity variables u, v, and w will be constant; however, this may not be trne 
for the density p and pressnre p. In the presence of gravity, a volnme of compressible 
finid will be compressed by the weight of the finid above it, increasing the density and 
pressnre. Therefore, it is reasonable to expect that onr reference states for density 
and pressnre will vary with 

Let ns consider a constant domain governed by (11.57). With everything at 
rest, we set all time derivative terms and all velocity terms to zero. This leaves ns 
with 


dxP = 0 
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dyP = 0 


(11.59) 


dzP = -gp 

This confirms our expectation that our reference values for p and p are dependent 
on So for our reference values, we let mq, vq, and wq define our constant mean 
velocities; p denotes our ^-dependent density reference state, with po the density at 
^ = 0 ; and p denotes our ^-dependent pressure reference state, with po the pressure 
at ^ = 0. 

2. Linearizing the Equations 

Begin with the hrst equation of (11.57): 

dtp + dx{pu) + dy{pv) + d:,{pw) = 0. (11.60) 

Expand the derivatives via the product rule, and substitute the reference/perturbation 
variable expansion in (11.58) to get 

(p + P*) + (p + P*) dx {uq + ) + {uq + u*) dx{p + p*) 

+ {p + p*) dy{vo + V*) + {vo + V*) dy{p + p*) >=0. (11.61) 

+ (p + p*) d^{wo + w*) + {wo + w*) d^{p + p*) 

/ 

Recalling which reference variables are independent in space and time, and eliminating 
all terms of 0(5^), we get 

dtp* + P {dxU* + dyV* + dz.w*) + u^dxP* + v^dyp* + w^d^p* = -d^p {wq + w*). (11.62) 

We leave the equation in this form, rather than reverting it to the whole state variable. 

By considering only the perturbation variable, we eliminate the possibility of the 
reference value overwhelming the perturbation, introducing unnecessary round-off 
errors into the hnite-precision calculations. 

When we apply this same process to the velocity and pressure equations of 
(11.57), using (11.59) to simplify the equation for w, we get 

dt<p + Adx(p + Bdy(p + Cdz<p = D(p + E, (11.63) 
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where 




A 


B 


C 


D 


E 


P* U* V* w* p* ^ 

^ uq p 0 0 0 ^ 

0 Mo 0 0 i 

0 0 Mo 0 0 

0 0 0 Mo 0 

^ 0 7P 0 0 Mo y 

^ Mo 0 p 0 0 ^ 

0 Mo 0 0 0 

0 0 Mo 0 4 

u p 

0 0 0 Mo 0 

^ 0 0 7P 0 Mo y 

^ tMo 0 0 p 0 ^ 

0 tMo 0 0 0 

0 0 tMo 0 0 

0 0 0 tMo 4 

^ 0 0 0 7P tMo y 

^ 0 0 0 -d,p 0 ^ 

0 0/00 
0-/0 0 0 

-c/4 0 0 0 0 

^ 0 0 0 0 y 

^ -dzpwo ^ 

fvo 

-fuQ 

0 

^ gPWo y 


(11.64) 
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For the ^{d^p, dyp, d^p} terms in the momentum equations, we use the following 
linearization: 


dxP 

P 


P 


dzP 

P 


dx {p + P* 
p + p* 
dxP* 


P \1 + ^ 


1 --+ - - 
P \P 


dxP* 

P 

dxP* 

P 

dyP* 

p 


d,p* + d,p 

p V p VP, 

dzP* PpA P*\ 


P 


P 


P 


dzP* I, P*\ 


where we use (11.59) in the next to last line. 

Note that the matrix A is singular if mq = 0 or Mq = 7 p/p. Similarly, B is 
singular if uq = 0 or Vq = ^p/p, and C is singular if wq = 0 or = 7 p/p. 

So what have we lost by this linearization? The main difference between non¬ 
linear flow and linear flow is the non-linear interaction between vortices [69]. For 
example. Fig. 4 shows a rising thermal bubble using non-linear equations (left) and 
their linearized form (right). The turbulence is clearly absent in the linearized case. 
However, wave motion is not noticeably affected. Fig. 5 shows an acoustic wave using 
the non-linear (left) and linearized equations (right). The differences are negligible. 
(These hgures were generated during some early work, applying a spectral element 
implementation of the non-linear system (11.55) in the xz plane influenced by gravity, 
and to the linearized form of the same system.) Since our non-reflecting boundary 
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conditions are intended for wave propagation, we do not need to keep the non-linear 
effects of vorticity. 




Fignre 4. A Rising Thermal Bnbble Using Non-linear (left) and Linear (right) Eqna- 
tions 




Fignre 5. An Aconstic Wave Using Non-linear (left) and Linear (right) Eqnations 

Before we begin developing NRBCs for this eqnation system, let ns spend some 
time discnssing NRBCs in general, with particnlar emphasis on the scalar-eqnation 
implementations of the NRBCs we are here adapting for a linear hrst-order system. 
This discnssion will be the snbject of the next chapter, and then the new NRBC 
implementations will be developed in Chapters IV-VI. 
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III. NON-REFLECTING BOUNDARY 

CONDITIONS 

A. OVERVIEW AND HISTORY 

We wish to solve fluid flow problems in only a portion of a large or infinite 
domain. By restricting our area of interest, we effectively create a boundary B where 
none exists physically, dividing our computational domain from the rest of the 
domain T). Implicit in our choice of B is the assumption that everything we wish to 
model is contained inside fl; nothing external impinges on our computational domain. 
The challenge we must overcome, then, is defining B in such a way that it behaves 
computationally as if there were no physical boundary. How, then, do we specify 
what occurs at these boundaries? The usual answer is to claim that waves flow out 
of the domain at the boundary, but they do not flow into the domain. If defined 
correctly, this claim will result in waves which propagate out of the computational 
domain without any reflection, so that the computational boundary is transparent 
to these outgoing waves, mimicking the real-world behavior where no such physical 
boundary exists. 

In general, there are two ways to simulate an open boundary. One may either 
surround the domain with an artificial absorbing medium, so that outgoing waves 
are diffused to zero before they return to the computational domain, or one may use 
a differential operator to prescribe the wave behavior at the boundary, so that only 
outgoing waves are permitted. 

Research into modeling open boundaries has been active since the late 1960s. 
Zienkiewicz and Newton [121] first derived the Sommerfeld radiation condition for 
outgoing wave propagation in 1969. In hindsight, this differential operator is surpris¬ 
ingly obvious and easy to derive. Beginning with the known general solution to the 
one-dimensional scalar wave equation 

dttu = c^dxxU , (III.l) 
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we have 


u{x, t) = F[x — ct) + G{x + ct) (III-2) 

for some functions F and G [105]. If F denotes the outgoing waves, and G the 
incoming, then we insist on G = 0 on the boundary. Differentiating u with respect 
to X and t, this gives us 

dtu = -cF' , d^u = F' . (III.3) 

From here we easily get the Sommerfeld boundary condition: 

dtU + cdxU = 0 . 

This boundary condition can be interpreted in two ways. The characteristic-based 
interpretation uses this condition as prescribing the Riemann invariant of the solution 
(see Ch. 8 of [24]). The wave-based interpretation describes it as requiring waves on 
the boundary to satisfy the one-way advection equation. Several early NRBCs were 
developed in the 1970s using these two interpretations. Wurtele et ah [119] used 
the characteristic method, while Pearson [93] and Orlanski [92] took the wave-based 
approach. Engquist and Majda [25, 26] extended the wave-based method, dehning 
a pseudo-differential operator solution to the 2-D wave equation and deriving Fade 
approximations thereto in a sequence of ever-more-accurate boundary conditions. 
Smith [101] took a simplistic, albeit computationally intensive, approach: Apply a 
Dirichlet boundary condition to one solution, apply a Neumann boundary condition 
to another, and then add the two solutions, cancelling the two reflections and leaving 
only the non-reflecting solution within the accuracy of the Neumann operator. One 
of the hrst absorbing layer methods was also published around this time by Davies 
[17] for the linearized Euler equations in a nested environment, using a “relaxation” 
function near the boundary to match the small-scale interior scheme with the large- 
scale outer model. Later NRBCs built on these methods or developed new approaches. 

In the 1980s, several new NRBC techniques appeared. Bayliss and Turkel 
[9, 10] developed a sequence of increasing-order boundary conditions based on an- 
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nihilating the first terms of the asymptotic expansion of the scalar wave eqnation. 
Miller and Thorpe [87] proposed a modihed Orlanski scheme based on alternative 
time-stepping methods. Perm and Gnstafsson [27] nsed Fonrier transformations to 
devise a downstream bonndary condition for the steady-state linearized Enler eqna- 
tions. Klemp and Dnrran [79] developed a “sponge layer” to absorb ontgoing gravity 
waves, applied to the linear Bonssinesq eqnations, later applied to the non-linear Enler 
eqnations by Giraldo and Restelli [32], Davies [18] compared varions techniqnes, in- 
clnding “diffnsion zones,” relaxation fnnctions, and radiation bonndary conditions. 
Raymond and Kno [95] modihed the Sommerfeld condition to consider tangential as 
well as normal derivatives in mnlti-dimensional hows. Ting and Miksis [110] pro¬ 
posed nsing Kirchhoh’s formnla to determine the bonndary valnes of waves exterior 
to a scatterer. Trefethen and Halpern [111] analyzed the Engqnist-Majda method; 
considering diherent approaches to approximating the psendo-diherential operator, 
they demonstrated its well-posedness and proved Engqnist’s and Majda’s proposition 
concerning which classes of Fade approximations are well-posed. Higdon [62, 64] de¬ 
veloped a seqnence of increasing-order bonndary conditions based on a prodnct of 
Sommerfeld conditions at varions incidence angles. Keller and Givoli [78] developed 
the “Dirichlet-to-Nenmann” (DtN) mapping, an NRBG method which converts the 
Dirichlet condition at inhnity to a Nenmann condition at the compntational domain 
bonndary; they then nsed this DtN mapping in a hnite element solntion of Laplace’s 
eqnation on an inhnite domain [38]. See also Givoli’s 1991 paper [33] reviewing the 
then-cnrrent state of the art. 

The 1990s saw an explosion in NRBG development. Kroner [82] adapted the 
Engqnist-Majda scheme to the 2-D linearized Enler eqnations, nsing Fonrier transfor¬ 
mations rather than psendo-diherential operators to dehne the bonndary condition; 
Giles [30] performed a similar techniqne, also nsing Fonrier transformations to de¬ 
rive an NRBG for the 2-D linearized Enler eqnations. G. Kreiss [80] nsed a simple 
Dirichlet condition at the downstream end for the pressnre term of the linearized 
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Euler equations, and he showed that the resulting error decreases exponentially up¬ 
stream of the open boundary. Collino [15] extended the Engquist-Majda scheme to 
open domains (requiring consideration of corner conditions) and to other wave-like 
equations. Tam and Webb [106] used the asymptotic expansion to dehne radiation 
boundary conditions compatible with their dispersion-relation-preserving hnite dif¬ 
ference scheme for the linearized Euler equations. Higdon continued his sequence 
of NRBC papers [65, 66, 67, 68], culminating in a robust NRBC sequence for the 
dispersive (Klein-Gordon) wave equation. Grote and Keller [47, 48, 49] extended the 
DtN technique to spherical waves and the Helmholtz equation, in hnite difference and 
hnite element methods; Thompson and Huan [109] later modihed this formulation to 
implement a hnite element solution of the spherical wave equation. Ren [96] used a 
2-D Sommerfeld-like boundary condition, 

{dt +c{axd^ +aydy))r] = 0 (HI.5) 

= 1 , 

to dehne an open boundary for the 2-D Boussinesq equations. Hagstrom and Hariha- 
ran [54] presented an NRBG for the 2-D and 3-D wave equation in polar/spherical co¬ 
ordinates, using auxiliary variables to remove the high-order normal derivative terms; 
Huan and Thompson [75] later modihed this scheme by using a spherical harmonic- 
based formulation. Safjan [98] also took the auxiliary variable approach, applying 
them to high-order Fade approximations to the pseudo-differential operator of the 
scalar wave equation. Jensen [77] compared several techniques, including NRBGs 
and sponge layers, for modeling open boundaries in a stratihed ocean model. 

The 1990s also saw the development of the Perfectly Matched Layer (PML), 
hrst dehned by Berenger [11] for the 2-D Maxwell equations. This absorbing layer 
method surrounds the computational domain with a dispersive medium, dehned in 
such a way that incoming waves at any incidence angle pass from the interior to 
the dispersive layer without any (theoretical) rehection. A 2007 paper by Skelton et 
al. [100] claimed to hnd over 1,000 references to Berenger’s paper in the literature. 
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This technique has been applied to the linearized Euler equations by Abarbanel et 
ah [1], Goodrich and Hagstrom [45], Hu [71, 72, 73], and Natal [88]; to the linearized 
shallow-water equations by Navon et ah [89]; to Maxwell’s equations using a second- 
order discretization scheme by Sjogreen and Petersson [99]; to the linearized and 
non-linear wave equation by Appelo and G. Kreiss [7] (following Appelo et ah’s well- 
posedness and stability theory in [6]); to elastic waveguides by Skelton et ah [100]; 
to the time-harmonic wave equation by Bermudez et ah [12]; and to the non-linear 
Euler and Navier-Stokes equations by Hu et ah [74], 

NRBG development continued in this decade as well, with new techniques and 
new applications of old techniques. Oliveira [91] combined a Sommerfeld condition 
with an absorbing layer to develop an NRBG for the transient mild-slope equation: 

V [CCgVr]) + {eCCg - a;2) 77 - dui] = 0 (HI.6) 

(see Eqn. (9) of [91]). Grote and Keller [50] developed an exact NRBG for the 3-D 
elastic wave equation with a spherical open boundary, based on annihilating the hrst 
terms of the wave’s spherical harmonics. Alpert et ah [4, 5] developed two NRBGs 
for the scalar wave equation, one using Hankel functions and Laplace transforms, and 
one using Fourier-Laplace transformations. Lie [83] used Fourier-Laplace transfor¬ 
mations to derive an NRBG for the shallow-water equations. Golonius and Ran [16] 
developed an absorbing buffer technique for conservation law-based systems by using 
Fourier transformations to filter the outgoing flow disturbances. McDonald [85, 86] 
derived a characteristic-based NRBG for the shallow-water equations and compared 
the performance of NRBGs and “relaxation zone” boundaries in nested-model envi¬ 
ronments. Hagstrom and Nordstrom [56] analyzed the use of extrapolation boundary 
values in solving the steady-state linearized Euler equations, showing the relation¬ 
ship between the position of the artihcial far-held boundary and the error norm of 
the discrete solution. Blayo and Debreu [13] considered a characteristic variable ap¬ 
proach to NRBGs in hrst-order systems for ocean and atmospheric models. Ghang 
et ah [14] used a Space-Time Gonservation Element and Solution Element (GE/SE) 
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method to solve the 1-D Euler equations, which effectively handled shockwaves inside 
the domain (although the results were good within the domain, the numerical results 
published showed poorer performance closer to the boundary). Guddati and him [51] 
used continued fractions (rather than the Engquist-Majda Fade approximations) to 
approximate the pseudo-differential operator, and they devised a formulation that 
could be applied to any convex polygon boundary (rather than the usual straight-line 
boundaries). Zahid and Guddati [120] incorporated PML-like “padding elements” 
into a continued fraction NRBG for dispersive waves. Atassi and Galan [8] used 
Fourier-Bessel modes to derive an NRBG for the non-linear Euler equations in an 
annular duct. Song and Bazyar [102] derived an NRBG for hnite element frequency- 
domain wave analysis based on Fade approximations. 

The 2000s also saw a revived interest in the Higdon NRBG sequence. Givoli 
and Neta created an algorithm to compute high-order hnite difference derivatives 
automatically, removing the algebraic complexity which limited the original Higdon 
sequence to third-order. This automation, along with an auxiliary variable method 
which removed the high-order normal derivatives, was applied to the Klein-Gordon 
equation and the shallow-water equations in a sequence of papers by them and their 
students and colleagues [39, 40, 41, 42, 43, 90, 113, 114, 115, 116]. Givoli again 
published a review of current NRBG techniques in [35]. In addition, Hagstrom and 
Warburton [57] developed a modihed form of the Givoli-Neta auxiliary variable NRBG 
for the scalar wave equation. This new method was expanded and analyzed in [37, 
53, 55, 84]. 

In this dissertation, we use the Givoli-Neta automation to apply the Hig¬ 
don NRBGs to the linearized Euler equations; we also extend the Givoli-Neta and 
Hagstrom-Warburton auxiliary variable NRBGs thereto. Some of the results pre¬ 
sented in the following chapters have been published or submitted for publication 
[19, 20, 21, 22]. The remainder of this chapter is devoted to the scalar form of the 
Higdon, Givoli-Neta, and Hagstrom-Warburton NRBG sequences, laying the ground- 
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work for their application to the linearized Enler eqnations. 


B. NRBCS FOR SCALAR EQUATIONS 
1. Higdon 

One of the simplest non-reflecting bonndary conditions (NRBCs) is the Som- 
merfeld radiation condition: 

{dt + c,d,)u = Q, (III.7) 

where u is the nnknown solntion to onr problem, and Cq is the wave propagation speed 
in the positive x direction (for a right-side bonndary; on other sides, replace with 
the appropriate normal derivative). In essence, this bonndary condition says that the 
ontgoing wave at the bonndary satishes the one-dimensional advection eqnation with 
advection speed Cq. This bonndary condition is most easily applied to the standard 
wave eqnation, 

duu = clV\ (III.8) 

It can also apply to the linearized shallow water eqnations (see [41]), 

dtu - fv = -gd^r] 

dtv + fu = -gdyTj (III.9) 

dt-q + ho {d^u + dyv) = 0 

whose snrface height component rj can be converted (see [113]) to the dispersive 
(Klein-Gordon) wave eqnation, 

dttV = clV'^rj - frj , co = \fgho (III. 10) 

The difficnlty in (III.9) and (III. 10) is that there is more than one wave speed. The 
same problem afflicts the standard wave eqnation in more than one dimension, as 
waves travelling in directions other than normal to the bonndary have wave speeds 
whose normal components are not eqnal to cq. Higdon proposed [68] dehning a 
bonndary condition as a prodnct of J Sommerfeld-like terms: 

j 

J\{dt + Cjd,)u = Q (III. 11) 
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Just as “a:?/ = 0” is a consolidation of “a: = 0 or ?/ = 0”, this product consolidates, into 
a single equation, numerous Sommerfeld conditions with different advection speeds. 
Since the differential operator is linear, it is easy to show that if any one of the 
individual Sommerfeld conditions is satished, then the consolidated Higdon condition 
is also satished. 

a. Reflection Coefficient 

Let us analyze the claim that this boundary condition is non-rehecting. 
Consider a wave-like equation with a solution of the form 

u(a;,t) (III. 12) 


This solution dehnes a wave with speed Cx- Applying the Higdon boundary condition 
to this equation, we claim that 

n = 0 (HI. 13) 

j=i 

If Cx does not exactly equal one of the Cj, this statement may not be true. To make 
it true, we must add a rehected wave, thus: 

j 

n + Cjdx) -f = 0 (HI. 14) 

i=i 

The magnitude of Rj is the reflection coefficient associated with the boundary con¬ 
dition. (If (HI. 13) is true, then Rj = 0.) Higdon [68] claims that this rehection 
coefficient is always less than one, and that it decreases as J increases. We present 
the proof here in more detail. 


Lemma III.l The magnitude of the reflection coefficient Rj of the Higdon NRBC 
defined in (III. 11), applied to a wave with speed Cx, is 


J 


= n 



C 


Proof. We prove by induction. For the J = 1 case, we have 


(HI.15) 


{dt J- cffix) 




0 , 


(HI.16) 
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and it is easy to show that 


|i?i| = ^ . (III.17) 

+ Cl 

Having proven there exists a J snch that this lemma is trne, we consider the J + 1 
case: 

J+i 

n (dt + c,d^) = 0 (111.18) 

i=i 

Since the differential operators are linear, we can write 

n (dt + C,d^) [(a* + = 0 . (III.19) 

i=i 

Expanding the derivatives inside the brackets and combining terms, we get 

n (dt + cjd,) = 0 . (III.20) 



(III.21) 


(III.22) 


Since these Cj and Cx are positive, each term in the prodnct is less than one. Thns, as J 
increases, the reffection coefficient gets smaller. Of conrse, well-chosen Cj’s will rednce 
the reffection coefficient more rapidly, bnt even poor choices (for example, failing to 
consider possible incidence angles or phase speeds) will still bring improvements. 
Hence, we can get good absorption either with a high J or with well-chosen c/s. 
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b. A Simplifying Assumption 

The difficulty in using (III. 11) comes from the rapidly increasing alge¬ 
braic complexity for large J. At J = 2, we have 




^2 


, ^ dtdx 


+ C2 


dtdx 


+ C1C2 


dx"^ 


u = 0 


(111.23) 


At J = 3, it becomes 


(4- r r r 

“I" "I" ^^dt^dx 


+ClC2g^ + ClCsg^ C2C3g^ + CiCsCs^) 


U 


y =0 


(III.24) 


and so forth. The Givoli-Neta algorithm automates the hnite difference computation 
of these high-order derivatives, but at a cost of 0{3'^) operations per time step. 
However, if we simply make all cj equal to a single wave speed c, then (III. 11) can be 
written as 

{dt + cd^yu = 0 (III.25) 


(see [90]), or, expanding the polynomial, 

f J 71 d'^ \ 

= 0 (111.26) 

We will show in the next section that this simplihcation, when applied to a hnite 
difference formulation, requires only O(J^) operations per time step. 

c. Discretization of Higdon NRBCs 

For our hnite diherence implementation of the NRBCs, we use hrst- 
order diherences for each derivative, with the diherence pointing into the domain and 
backward in time. Higdon [68] demonstrated that this discretization scheme is stable 
when used in conjunction with the Klein-Gordon equation. For an open boundary on 
a rectangular domain, we have 


30 



0 


(111.27) 


Bottom: o", = 0, 

where I denotes the identity operator; S~ and S’’*' backward and forward shifts in the 
snbscript variable, e.g., 


Q+ n _ n 


Q—^77 _ n—1. 

^i,j — 5 


St onr grid spacing in time; Sx and Sy onr grid spacing in x and y, respectively; afj 
the valne of a generic state variable a at grid point {i,j) at time step n; and snbscripts 
N, W, E, S denoting the north, west, east, or sonth bonndaries, respectively. If we 
mnltiply each term by St, clearing it ont of the denominator, and then gronp terms 
by shift operator, we get 


where 


Top: (^ttyl + bS^ + CySy^ = 0 

Left: (aJ + bS;~ + CxS^^ = 0 

Right: (a^I + bS^ + CxS~^ = 0 

Bottom: (ayl + bS^ + CySy ^ a'^g = 0, 


(III.28) 


ttr = 1 — Cr 


Qjy 1 Cy 


b = 


('X Cq 


Cy 


a 

Sx 
St 


Expanding these operators as polynomials, we have 

Top: IE E 


cr. 


i,N 


0 
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Left: 

Right: 

Bottom: 


j J-0 


EE 

3=0 7=0 
( J J-0 

EE 


J! 




J! 


J J-0 

EE 


a. 




E,j 


J! 


fr«c2S,-'’sr < 


0 

0 


(111.29) 


where a = J — (3 — 'y. 

Note. From this point forward in the dissertation, we will only show 
the NRBC formnla for one side rather than all fonr. On the other three sides, the 
appropriate changes shonld be made for the correct normal derivative. 


2. Givoli-Neta 

The greatest problem afflicting these high-order NRBCs is the presence of 
high-order spatial and temporal derivatives. With the spatial derivatives, the NRBC 
algorithm mnst look deep into the domain, and an incoming wave can begin to affect 
the bonndary long before it actnally reaches the bonndary. High temporal derivatives 
reqnire a long “history” of past valnes, increasing memory reqnirements. In both 
cases, and with the high-order mixed derivatives, increasing the nnmber of terms in 
the NRBC calcnlation increases the danger of ronnd-off errors corrnpting the solntion 
and destabilizing the system. 

In addition, the NRBC order is inherently limited by the size of the domain. 
This is trne in a hnite difference setting, bnt it is even more restrictive in dealing with 
hnite elements. In a finite element scheme, the spatial derivatives generally limit the 
NRBC order to that of the element polynomials, u nl ess one is willing to nse derivative 
approximations which are not local to a single element. 

The way ont of this dilemma is to introdnce auxiliary variables [34] which 
reqnire only low-order derivative calcnlations. In this section and the one following, 
we consider two methods for implementing these anxiliary variable techniqnes. 

The Givoli-Neta NRBC was hrst proposed in [39, 42, 43] for the Klein-Gordon 


32 



equation in finite difference and finite element schemes. (In this dissertation, the 
term “Givoli-Neta NRBC” refers to the auxiliary variable implementation, not their 
automated high-order Higdon NRBC.) This NRBC follows directly from the Higdon 
NRBC. Where the Higdon NRBC uses one high-order equation on the boundary, the 
Givoli-Neta NRBC of order J uses a system of J low-order equations, dehned thus: 

^j+i = ‘Pi (HI.30) 

(po = u 

(pj = 0 . 

Direct substitution shows that (HI.30) is equivalent to (HI. 11). The above dehnition 
applies to an NRBC on the right side of the domain; on other sides, replace dx with 
the appropriate normal derivative. 

To illustrate the utility of this formulation, let us apply the NRBC to the 
Klein-Gordon equation (HI. 10). As this work has already been done by Givoli and 
Neta [39, 42, 43], we will not repeat the derivations here. 

It is easy to show that the auxiliary variables also satisfy the Klein-Gordon 
equation, that is, 

dtt^j = CoVVi - /Vi , (HI.31) 

for all J G 1... J — 1. Armed with that fact, we can then combine (HI.30) and (HI. 10) 
to remove the normal derivative terms from our NRBC. When we do so, our resulting 
NRBC is 

( 72 ~ 3 ) A [ “ I" ~ I dtp>j — p>j+i = dyytpj-i — -^(pj-i . (HI.32) 

\Co C/ VC C+i/ C 

This equation gives us a tri-diagonal system of equations to solve for our 
auxiliary variables. With an efficient matrix solver, the effort required is 0{J), rather 
than the O(J^) required for the original formulation, and it is 0(J) even without 
setting the Cj values to Cq. Furthermore, since the auxiliary variables are dehned solely 
using temporal and tangential derivatives, we can restrict them to the boundary. 


33 



reducing the amount of memory needed to store their values. Finally, since the 
NRBC uses no derivatives beyond second order, they can be easily applied to a hnite 
element system to any order. There are three downsides, however. NRBCs on two 
adjacent sides require consideration of corner compatibility conditions (due to the 
presence of tangential derivatives), their numeric stability is less robust over long time- 
integrations [37], and the conversion to remove the normal derivatives introduces an 
additional source of discretization error which produces an “error floor” that cannot 
be overcome by simply increasing J [36]. 

3. Hagstrom-Warburton 

The Hagstrom-Warburton NRBCs were hrst presented in [57] for the scalar 
wave equation, with subsequent extensions and analysis in [37, 53, 55]. This NRBC 
scheme also uses a system of low-order auxiliary variable equations instead of a single 
high-order boundary equation. The Hagstrom-Warburton NRBC of order J is given 
by 


dt^pi = aodtu + coda;U 
ajdt(pj+i - cod^ipj+i = Qjdtipj + cod^Lpj 
^j+i = 0 , 


(HI.33) 


where the parameters aj G (0,1] are chosen by the user. Hagstrom and Warburton 
show that this dehnition, applied to a plane wave traveling at speed cq and incidence 
angle 9, results in a reflection coefficient of 

Oo — cos 9 -^ ( aj — cos 6* \ ^ 
oo + cos 9 y aj -T cos 9 j ’ 


R = 


(HI.34) 


which, like the Higdon scheme’s reflection coefficient, is a product of J terms, each 
of which is less than unity. Hence, as J increases, the reflection coefficient decreases. 
However, with the squaring of each term in the product, the reffection coefficient 
for the Hagstrom-Warburton scheme decreases signihcantly more rapidly than the 
Higdon reflection coefficient. 
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No physical interpretation of the NRBC formulation is given. One way to 
interpret (III.33b) is that the outgoing characteristic of (pj is matched by the incoming 
characteristic of pj+i. However, that does not explain the form of (III.33a). It could 
be that the authors wanted an NRBC scheme with a quadratically-decaying reflection 
coefficient, and they reverse-engineered this NRBC formulation to get one that works. 

As with the Givoli-Neta NRBCs, the auxiliary variables presented here also 
satisfy the Klein-Gordon equation (III.31) for all j G 1... J. We can use that fact to 
remove the normal derivative term from (III.33b), replacing that equation with 


- l) dtt<fij-i + aj-i (aj - l) duPj+i 
— (aj_i -f- ttj) {aj-iaj -h 1) duPj 


> = < 


(/Vi-l - cldyyPj-l) 

+ (%-l + %) (/Vj - Cldyypj) 

+aj_i (/Vi+1 - cldyyPj+i) 

(III.35) 

(again, see [55] for details). Once again, we have a system of auxiliary variables 
dehned solely on the boundary, and a boundary condition consisting entirely of low- 
order derivatives. As with the Givoli-Neta system, this system also requires 0{J) 
operations to solve. 


C. NRBCS FOR FIRST-ORDER SYSTEMS 

In general, the three NRBC techniques discussed in the previous section have 
only been implemented for scalar wave-like equations, not for hrst-order systems. 
The exception is the Hagstrom-Warburton NRBC, where one section of [57] briefly 
discusses the characteristic-based implementation for a symmetric hrst-order system 
of the form 

dtp + AdxP + Bdyp = 0 . (III.36) 

However, that discussion does not consider corner conditions or the presence of un¬ 
differentiated terms on the right-hand side. In Chapter VI we develop these NRBCs 
in more depth. 
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While the Higdon and Givoli-Neta NRBC methods have been developed for 
the shallow-water eqnations (e.g., [113]), those implementations involve converting 
the system to a scalar Klein-Gordon eqnation for the snrface height state variable; 
hence, they are not trnly implemented for the hrst-order system. In the next two 
chapters, we will develop these NRBG methods for a trne hrst-order system. While 
the Higdon implementation is specihc to the linearized Enler eqnations, the Givoli- 
Neta implementation (not inclnding gravitational effects (Sec. V.G)) can be extended 
easily to any hrst-order system. 
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IV. HIGDON NRBCS FOR THE 
LINEARIZED EULER EQUATIONS 


A. OVERVIEW 

This chapter develops the hnite difference implementation of the high-order 
NRBCs for the linearized Enler eqnations in two dimensions. We begin with the 
simplest possible eqnation set, with no advection or forcing terms. After demon¬ 
strating the validity of this prototypical implementation in Section B, we proceed to 
incorporate the effects of Coriolis (Sec. C), gravity (Sec. D), and advection (Sec. E). 


B. AN INITIAL PROTOTYPE 

We begin with the simplest possible implementation of the linearized 2-D Enler 
eqnations: zero advection, no Coriolis or gravitational forces. 


dtP + Po + dyv) = 0 

dtu + —dxP = 0 
Po 

dtv + —dyp = 0 

Po 

dtp + 7Po {dxU + dyv) = 0 


(IV.l) 


The variables here are the pertnrbation variables from the linearization in Sec. II.C.2; 
we remove the asterisks here and in all snbseqnent eqnations for textnal clarity. In 
addition, we nse po and po instead of p and p becanse we have no ^-dependency in these 
eqnations. This set is mostly nseless for real-world modeling, especially considering 
the fact that p has no impact on any of the other state variables. However, it will 
serve as an initial prototype for developing this implementation (see also [19]). 

1. Equivalence of (IV.l) and the Scalar Wave Equation— 
Continuous Case 

In [68], Higdon proved that the J-order NRBC is compatible and stable when 
applied to the dispersive wave eqnation (and thns the standard wave eqnation by 
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setting / = 0), based on stability criteria developed by H.-O. Kreiss [81] (with a 
geometric characteristic-based interpretation provided by Higdon [63]). Therefore, if 
we can convert these simplihed Enler eqnations to the standard wave eqnation, we 
will know that they too will be stable with this NRBC formnlation. 

Differentiate (IV.lb,c,d) with respect to x, y, and t, respectively, and get 
(ignoring the eqnation for p) 

dxt'U' dxxP 

Po 

dytV = -—dyyP (lV.2) 

Po 

dttP = -IPO {dxtu + dytv) 


Snbstitnting the hrst and second of these new eqnations into the third gives 

dttP = — {d^xP + dyyp ), 

Po 

which is the standard wave eqnation for p, with the wave speed given by 


(IV.3) 


Co = 


npo 

Po ’ 


(IV.4) 


reflecting the valne given in the literatnre [24, 76, 107, 108]. Therefore, this system of 
eqnations, combined with the Higdon-like NRBCs, will be stable. It may seem a bit 
shaky to make this claim, since we have only proven that the system can be converted 
to the wave eqnation for p, not for the other variables. Might they still be nnstable 
in u and n? The answer is no. Since p depends on u and n, if they are nnstable, 
then they will destabilize p. Since p is stable, it follows that u and v mnst also be 
stable. (The stability of p is somewhat irrelevant, as it has no impact on the other 
state variables.) 

In Sec. VI. 1, we will show that all fonr state variables satisfy the scalar wave 
eqnation, and we will derive its exact form. For now, it is snfficient to note that we 
can convert this eqnation set to the scalar wave eqnation in p. 
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2. Equivalence of (IV. 1) and the Scalar Wave Equation— 
Discrete Case 

Using stability criteria developed by Gustafsson et.al. [52], Higdon showed 
[68] that the discretization scheme in Sec. III.B.l.c is compatible with the standard 
second-order discretization scheme for the standard wave eqnation (solving for 


u. 


n+l 


2«m- 




u. 


{Sty 


= Cn 




2«m- 


+ <-i, 


(fe) 


+ 


u 




2<. 


+Ki-y 


{Syy 


(IV.5) 


However, this scheme is for solving a single second-order PDE, not a system of hrst- 
order PDEs, as we have. Even thongh we have shown that this system is eqnivalent, 
in one of the state variables, to the standard wave eqnation, creating a discretization 
scheme which is also consistent and stable is a matter of some delicacy. 

Givoli and Neta described a similar difficnlty they enconntered in developing 
a stable discretization scheme for the shallow water eqnations [39]. After nnmerons 
failed efforts, whose resnlts were nnstable for any J > 2, they dehned a scheme which 
was eqnivalent to the dispersive wave eqnation discretization scheme proved stable 
by Higdon in [68]. We nse a similar approach here. Let denote a second-order 
centered difference in time, with similar notation for difference approximations in x 
and y; thns. 


A 


a 


a 


s* - s; 

25a 

e {x,y,t} , 


(IV.6) 


where 5a is the step size in a. Use the following discretization for (IV.1): 


Atp 


Atti 


/\tv 




-po {A^u + Ayv) 
A^p 


Po 

\P 

Po 

-7Po {A^u + Ayv) 


(IV.7) 
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As with converting (IV. 1) to (HI.8), apply A^. to (IV.7b), to (IV.7c), A^ to (IV.7d), 

and then make the appropriate snbstitntion. This gives ns 



which is the scalar wave eqnation on a double-sized grid. Hence, to make onr dis¬ 
cretization scheme fnlly compatible, we nse the following discretization for the NRBC 
on each side of the domain: 






(IV. 10) 


where a denotes any one of onr state variables p, u, v, p. 

A reviewer for [21] noted that this NRBC cannot resolve the shortest wave¬ 
lengths resolvable by the interior scheme. Snbseqnent experiments showed no new 
instabilities generated by these short wavelengths. 

In [19], the anthors showed that a certain choice of one-sided differencing can 
also yield a compatible scheme, one which is compatible with the scalar wave eqna¬ 
tion on the same size grid. However, snbseqnent experimentation revealed deeper 
flaws in that approach. First, it conld not be extended to match the Klein-Gordon 
eqnation when Coriolis forces are inclnded. Second, its semi-implicit natnre became 
fnlly implicit nnder advection. (An attempt to avoid an implicit scheme, based on 
a qnadrant-based interior discretization method, proved to be nnstable nnder ad¬ 
vection.) The scheme devised here avoids these flaws, as snbseqnent sections will 
show. 


3. Numerical Example—Semi-Infinite Channel 

Consider a simple nnmerical example. Dehne a sqnare domain 10 km on each 
side, with walls on three sides and open on the top (see Fig. 6). We trnncate the 
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domain with an artificial north boundary T and impose the non-reflecting boundary 
condition on the state variables. We dehne the hard wall as 


dnP = 0 

u-n = 0 (IV.ll) 

dnP = 0 

The open boundary is dehned by a Higdon-type NRBC of order J. Using a sea-level 



Figure 6. A semi-inhnite channel domain truncated by an artihcal boundary F^ 
[After [40], Fig. lb, p. 259] 

atmospheric density po = 1-2^ and pressure po = 1.01 x 10®^ [58], our initial 
condition is a pressure bubble in the center of the domain: 


Px,z = 


To 1 + 


100 


d < r 


Po : otherwise 


pI,z 


= < 


Po 


4^1" : d<r 


Pz 


Po : otherwise. 


(IV.12) 
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where 



and {xc, Vc) denotes the center of the domain. The initial pertnrbation of p is dehned 
to maintain a constant potential temperatnre in the pressnre pertnrbation bnbble 
[24], We divide the domain into a 100 x 100 grid and rnn the simnlation np to t = 24 
s, which is snfficient time for the wave trongh to reach the corners. We also rnn 
a reference solntion on a 10 x 20 km domain consisting of 100 x 200 grid points; 
this reference solntion has hard walls on all fonr sides, and the simnlation dnration 
is short enongh that reflections from the top bonndary will not re-enter the original 
domain. By nsing a compnted reference solntion rather than an analytic solntion, we 
can attribnte all differences to the NRBC error rather than the discretization scheme’s 
trnncation error. Onr time step is chosen as 90% of the CFL limit [64], 




2 


< 1 , 


(IV.13) 


where Cq is the wave speed \J'^Po/po- (We have observed experimentally [21] that 
setting 5t to the maximnm allowed by the CFL limit resnlts in rednced effectiveness 
for the higher J NRBCs.) We nse this same cq for the NRBC wave speed. Dehne the 
error norm, for each state variable p, to be 


E^p — 




(IV. 14) 


where pj is the state variable compnted nsing the order J NRBC; pq is the reference 


solntion state variable; W and Ny are the nnmber of points in the x and y directions, 
respectively; and we normalize the error norm by the norm of the reference state 
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variable (so that all four state variables’ error norms are approximately the same 
order of magnitude). Using the perturbation variables dehned in (IV. 1), the interior 
discretization scheme (IV.7), and the Higdon NRBC discretization (IV. 10a), we run 
the simulation using different values of J for the NRBC order. Figs. 7-10 show 
the perturbation variables p, u, u, and p, respectively, for the J = 10 case. Fig. 11 
contrasts the computed solutions and error deltas for v between the J = 1 and J = 10 
cases. Table I shows the error norms for each state variable for J G 1... 10. 



Solution delta 

100 
do 
60 
40 
20 

50 100 

Error ftorm-0.0022101 



A 

X 10 



Figure 7. Plot of p in basic system (IV. 1) with J = 10 in a semi-inhnite channel. 
(TL) Computed solution. (Center) Reference solution; the area corresponding to the 
computed solution is contained below the horizontal line. (BL) Reference solution 
truncated to computed solution domain. (BR) Delta between reference solution and 
computed solution, with error norm computed by (IV. 14). 
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Figure 8. Plot of u in basic system (IV. 1) with J = 10 in a semi-infinite channel. 
(TL) Computed solution. (Center) Reference solution; the area corresponding to the 
computed solution is contained below the horizontal line. (BL) Reference solution 
truncated to computed solution domain. (BR) Delta between reference solution and 
computed solution, with error norm computed by (IV. 14). 


J 

Ep 

Eu 

Ey 

Ep 

1 

0.12361 

0.077449 

0.1674 

0.12361 

2 

0.039728 

0.020397 

0.047278 

0.039727 

3 

0.026079 

0.010879 

0.022051 

0.026078 

4 

0.022763 

0.0077222 

0.014192 

0.022763 

5 

0.021505 

0.0060029 

0.010178 

0.021505 

6 

0.020889 

0.0049066 

0.0078822 

0.020889 

7 

0.020551 

0.0041032 

0.0063328 

0.020551 

8 

0.020362 

0.0035019 

0.0052668 

0.020361 

9 

0.020249 

0.0030695 

0.0044665 

0.020249 

10 

0.020176 

0.0027737 

0.0038836 

0.020176 


Table I. Error norms for basic system (IV.1) with J G 1... 10 in a semi-inhnite 
channel 
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Figure 9. Plot of v in basic system (IV. 1) with J = 10 in a semi-infinite channel. 
(TL) Computed solution. (Center) Reference solution; the area corresponding to the 
computed solution is contained below the horizontal line. (BL) Reference solution 
truncated to computed solution domain. (BR) Delta between reference solution and 
computed solution, with error norm computed by (IV. 14). 
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Figure 10. Plot of p in basic system (IV.1) with J = 10 in a semi-infinite channel. 
(TL) Computed solution. (Center) Reference solution; the area corresponding to the 
computed solution is contained below the horizontal line. (BL) Reference solution 
truncated to computed solution domain. (BR) Delta between reference solution and 
computed solution, with error norm computed by (IV. 14). 
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Figure 11. Comparison of v in basic system (IV. 1) computed with J = 1 and J = 10 
in a semi-inhnite channel, with error norms computed by (IV. 14). (TL) Computed 
solution for J = 1. (TR) Computed solution for J = 10. (BL) Delta between reference 
solution and J = 1 computed solution. (BR) Delta between reference solution and 
J = 10 computed solution. 
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4. Numerical Example—Open Domain 

We now perform the same simulation in an open domain (see Fig. 12), using 
NRBCs on all four sides for the state variables. Before proceeding, we note that we 


Figure 12. An open domain 12 truncated by artihcial boundaries Fl, F^, Fi^ and F b 
[After [19], Fig. 1, p. 1] 

now have points which are on the corners of two open boundaries, and we pause to 
address this potential complication. Looking at (IV. 10), we see that for each point 
on an open boundary, that point only depends on itself at previous times and on 
the adjacent interior points. The boundary points are independent of each other. 
Furthermore, due to the discretization (IV.7) of the interior scheme, there are no 
interior points which depend on the corner points. Fig. 13 illustrates this dependency 
for the top-right corner of an open domain. Notice that no other point depends on 
the corner. Hence, our boundary condition at the corner is, by and large, irrelevant. 
For our implementation here, we decree that the corners are considered part of the 
top/bottom boundaries. Our reference solution is a 30 x 30 km domain with 300x300 
grid points, situated so that the center of the reference domain corresponds to the 
computational domain. When we run this simulation, we get results such as those in 
Figs. 14, 15, and Table II. (From here on, we will only plot the results for one state 
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Figure 13. Interior and boundary discretization dependencies for points near the 
top-right corner of an open domain. The black arrows show the interior dependecies 
based on the discretization scheme (IV. 7). Blue arrows show the dependencies of the 
boundary points except for the corner. Green arrows show the dependency if the 
corner is considered part of the top; red arrows show the dependency if the corner is 
considered part of the right. 

variable, rather than all four.) Not surprisingly, given the symmetry of the domain 
and the discretization scheme, the errors for u and v are the same. 


J 

E, 

Eu 


Ep 

1 

1.5544 

2.0918 

2.0918 

1.5558 

2 

0.4589 

0.61777 

0.61777 

0.45933 

3 

0.22861 

0.30055 

0.30055 

0.22882 

4 

0.15198 

0.19766 

0.19766 

0.15212 

5 

0.11326 

0.14588 

0.14588 

0.11336 

6 

0.089796 

0.11564 

0.11564 

0.08988 

7 

0.074067 

0.095798 

0.095798 

0.074136 

8 

0.06246 

0.082285 

0.082285 

0.062519 

9 

0.053427 

0.071617 

0.071617 

0.053476 

10 

0.046644 

0.062478 

0.062476 

0.046687 


Table II. Error norms for basic system (IV. 1) with J G 1... 10 in an open domain 
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Figure 14. Plot of u in basic system (IV.1) with J = 10 in an open domain. (TL) 
Computed solution. (Right) Reference solution; the area corresponding to the com¬ 
puted solution is contained within the center box. (CL) Reference solution truncated 
to computed solution domain. (BL) Delta between reference solution and computed 
solution, with error norm computed by (IV. 14). 



Figure 15. Comparison ofp in basic system (IV.1) computed with J = 1 and J = 10 
in an open domain, with error norms computed by (IV. 14). (TL) Computed solution 
for J = 1. (TR) Computed solution for J = 10. (BL) Delta between reference 
solution and J = 1 computed solution. (BR) Delta between reference solution and 
J = 10 computed solution. 
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C. CORIOLIS FORCES IN THE XV PLANE 

Let us now add Coriolis forces to our equation set [21], Leaving the advection 
terms zero, our equation set is 


dtP + Po (Xu + dyv) = 0 
dtu + —dxP = fv 


Po 

dtv H- dyp = -fu 

Po 

dtp + '^Po(d^u + dyv) = 0 (IV. 15) 

This set is somewhat more useful for atmospheric modeling, although we still have p 
decoupled from the rest of the system. As with the preceding section, we now wish 
to show that (IV. 15) is equivalent to the Klein-Gordon equation duP = CqV^P — /^p, 
which, as noted previously, Higdon proved is stable with this NRBC formulation [68]. 


1. Equivalence of (IV. 15) and the Klein-Gordon Equation— 

Continuous Case 

This conversion begins the same as in the previous section. Differentiate 
(IV.15d) with respect to t 


dttP + 'fPoiXtU + dytv) = 0 


(IV.16) 


Now differentiate (IV. 15b) with respect to x and (IV. 15c) with respect to y and add 


^xtU T OytV (dxxP T ^yyP^ f (^xU OyU^ 
Po 


(IV. 17) 


Now substitute (IV. 17) into (IV.16) 


dttP - — V^p + /7Po (da^v - dyu) = 0 . 
Po 


(IV.18) 


Differentiate (IV. 15b) with respect to y and (IV. 15c) with respect to x and subtract 

( 


OytU OxtU T 

Po 


9yxP dxyP f (dyV T dxU^ 


(IV.19) 


V 


=0 
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Combine terms to get 


/ {d^u + dyv) = -dt {d^v - dyu) . (IV.20) 

Combine (IV.15d) and (IV.20) to get 

fdtp = 'ypodt {da;V - dyu) . (IV.21) 

Integrate (IV.21) with respect to time to get 

fip- Po) = IPo {dxV - dyu) . (IV.22) 

Finally, substitute (IV.22) into (IV. 18) 

dttP - —V^p + f{p-po) = 0, (IV.23) 

Po 

which gives us the Klein-Gordon equation for the pressure perturbation p — po, again 
with wave speed \j'ypo/Po- As an aside, we notice that (IV.20) involves the time 
derivative of the curl of u on the right-hand side. Thus, if / = 0, then V x F is 
constant. 


2. Equivalence of (IV. 15) and the Klein-Gordon Equation- 
Discrete Case 

We would like to hnd a discretization scheme which is equivalent to the discrete 


Klein-Gordon equation 


_ 27/” -I- 7;”t 

J 


.n—1 


U. 


St'^ 


= Cn 


r+i.i - 2<i+<-ij , <i+i “ + Kj-i 


SI 


+ 
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(IV.24) 

Given the similarity to the discrete scalar wave equation, we begin by using the 
discretization scheme in Sec. B.2, adding the Goriolis terms to the right-hand side to 
get 

^tP = -po (AxM + ^yV) 

Atu = -^^ + fv 


Po 

\P 

Po 


fu 


= -7Po (A^u + Ayv) 


Atv 

Atp 
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(IV.25) 



Apply Aj. to (IV.25b), to (IV.25c), A^ to (IV.25d), and then make the appropriate 

snbstitntion. This gives ns 


AtAtp = — {A^A^p + AyAyp) - f'fpo {A^v - Ayu ), (IV.26) 

Po 

Next, apply Ay to (IV.25b) and A^, to (IV.25c), then snbtract and combine terms to 
get 

At {AyU - A^v) = f {A^u + Ayv) (IV.27) 

We snbstitnte (IV.25d) into (IV.27) to get 

At {AyU - A^v) = -^Atp (IV.28) 

If we apply At to (IV.26) and incorporate (IV.28) into it, we get 


At AtAtp - {A^A^p + AyAyp) + f p =0 

Po 


(IV.29) 


Thns, the qnantity inside the brackets is constant in time. Since we are assnming 
a closed system with no sonrce fnnctions (other than the Coriolis force itself), then 
this qnantity mnst initially be zero and will thns always be zero. Expanding the 
differencing terms gives ns 

^ in ( p^2,-^p?,+pi2, , , 

(25t)^ Po \ {2SxY (25|/)^ / 

(IV.30) 

As with eqnating the initial basic case to the scalar wave eqnation, we now have a 
system which is compatible with the Klein-Gordon eqnation on a donble-size grid. So 
again, onr NRBC will be given by (IV. 10). 


3. Numerical Examples 

For this example, we nse the same domain and initial conditions as with the 
basic system in Sec. B. Onr Coriolis force / is based on an angnlar momemtnm 
O = 7.292116 X 10“® s“^ [112] at a latitnde (j) of 30° N, thns. 


/ = 20sin()) = 7.292116 x 10“® s“^ 
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J 

E, 

Eu 

Ey 

Ep 

1 

0.12361 

0.077448 

0.1674 

0.12361 

2 

0.039728 

0.020397 

0.047278 

0.039727 

3 

0.026079 

0.010879 

0.022051 

0.026078 

4 

0.022763 

0.0077222 

0.014192 

0.022763 

5 

0.021505 

0.0060029 

0.010178 

0.021505 

6 

0.020889 

0.0049066 

0.0078821 

0.020889 

7 

0.020551 

0.0041032 

0.0063328 

0.020551 

8 

0.020362 

0.0035019 

0.0052668 

0.020361 

9 

0.020249 

0.0030695 

0.0044665 

0.020249 

10 

0.020176 

0.0027737 

0.0038836 

0.020176 


Table III. Error norms for Coriolis-influenced system (IV. 15) with J G 1... 10 in a 
semi-infinite channel 


When we apply this force to our domains, we get results which are virtually indistin¬ 
guishable from those in the preceding section. This is not too surprising, since the 
Coriolis force is so small that its effects are unnoticeable over such short timeframes 
[28]. Tables III and IV show the error norms for all four state variables for J G 1... 10 
in the channel and open domain, respectively. Suppose, just for curiosity’s sake, we 
make the Earth rotate much more rapidly and increase / by a factor of a thousand. 


J 

E, 

Eu 

Eu 

Ep 

1 

1.5544 

2.0917 

2.0917 

1.5558 

2 

0.4589 

0.61777 

0.61777 

0.45933 

3 

0.22861 

0.30055 

0.30054 

0.22882 

4 

0.15198 

0.19766 

0.19766 

0.15212 

5 

0.11326 

0.14588 

0.14588 

0.11336 

6 

0.089797 

0.11564 

0.11564 

0.08988 

7 

0.074067 

0.095798 

0.095797 

0.074136 

8 

0.062461 

0.082285 

0.082284 

0.062519 

9 

0.053427 

0.071617 

0.071617 

0.053476 

10 

0.046645 

0.062478 

0.062476 

0.046689 


Table IV. Error norms for Coriolis-influenced system (IV. 15) with J G 1... 10 in an 
open domain 
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J 

E, 

Eu 

Ey 

Ep 

1 

1.4093 

0.70281 

0.70253 

1.4033 

2 

0.4086 

0.20314 

0.20309 

0.40687 

3 

0.20321 

0.098226 

0.098178 

0.20235 

4 

0.13492 

0.064387 

0.064365 

0.13435 

5 

0.10056 

0.047464 

0.047443 

0.10013 

6 

0.079742 

0.037609 

0.037594 

0.079405 

7 

0.065817 

0.031154 

0.031144 

0.065539 

8 

0.055543 

0.026777 

0.026767 

0.055308 

9 

0.047521 

0.023331 

0.023321 

0.04732 

10 

0.041445 

0.020369 

0.020361 

0.04127 


Table V. Error norms with J G 1... 10 in an open domain nnder artifically-large 
Coriolis (IV. 15) 


to / = 0.07292116 s“^, jnst to make its effects noticeable. If we do that, we get the 
error norms given in Table V for the open domain. A comparison of the resnlts for 
small and large J is exemplified in Figs. 16 and 17. If we look at the u and v plots 
side-by-side (Fig. 18), we can discern the clockwise rotation generated by the Coriolis 
force. 
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X-velocity perturbation, J^l 



X-vek>city perturbation error, Jsi 



I 

I 

I 

0 

In -( 


X*velocity perturbation. JslO 





0.1 

0.05 

0 

-0.05 

• 0,1 


X-velodty perturbation error. J^IO 



I 

0 
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Figure 16. Comparison of u in Coriolis-influenced system (IV. 15), using artificially- 
large Coriolis, computed with J = 1 and J = 10 in an open domain, with error norms 
computed by (IV. 14). (TL) Computed solution for J = 1. (TR) Computed solution 
for J = 10. (BL) Delta between reference solution and J = 1 computed solution. 
(BR) Delta between reference solution and J = 10 computed solution. 


Y-velocity perturbation, J=1 Y-velocity perturbation. J=10 



20 40 60 80 100 20 40 60 80 100 



Figure 17. Comparison of v in Coriolis-influenced system (IV. 15), using artihcially- 
large Coriolis, computed with J = 1 and J = 10 in an open domain, with error norms 
computed by (IV. 14). (TL) Computed solution for J = 1. (TR) Computed solution 
for J = 10. (BL) Delta between reference solution and J = 1 computed solution. 
(BR) Delta between reference solution and J = 10 computed solution. 
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Figure 18. Side-by-side plot of u and u, illustrating the rotation generated by the 
Coriolis acceleration. The superimposed clockwise arrows highlight the combined 
result of the u and v values. 
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D. GRAVITATIONAL FORCES IN THE AZ PLANE 


Now we shift from the xy plane to the xz plane and consider what happens 
when we add the effects of gravity to onr system. Still keeping the advection terms 
eqnal to zero, we now have the eqnation set 


dtp + p {d^u + d^w) = -p'w 

dtu + -d^p = 0 

P 

dtw + -d^p = --p 

P P 

dtP + 7P (Am + d^w) = gpw 


(IV.31) 


where we now mnst nse p and p for onr ^-dependent reference states. We nse p' to 
denote the derivative of p, which depends only on Note that p hnally makes an 
impact on the other state variables. 

1. Defining the Reference State for Density and Pres¬ 
sure 

Eq. (11.59) sets a compatibility condition between onr reference states for den¬ 
sity and pressnre. This restriction is a fairly loose one, however, reqniring ns to look 
to other sonrces for possible fnnctions which satisfy this condition. Since these eqna- 
tions are derived from the ideal gas law, we look to the literatnre for atmospheric 
models on which we can base these initial conditions. 

Althongh several atmospheric models exist [24, 112], for ease of differentiation 
we will nse an exponentially-decaying model 


— — OlZ 

p = Poe 


(IV.32) 


where a is a scaling height needed to match the snrface {z = 0) pressnre and density 
valnes. Applying (11.59) to this model, we get 


P 


Pqo: 

-e 

9 


(IV.33) 
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2. Deriving a Wave-Like Equation from (IV.31) 

As with the basic system and the system subject to Coriolis forces, we now 
attempt to convert the gravity-affected system to a wave-like equation. First, we 
differentiate (IV.31d) with respect to time to get 

dttP = gpdtw - 7p {d^tu + d^tw) (IV.34) 


Differentiating (IV.31b) with respect to x and (IV.31c) with respect to ^ (remember¬ 
ing that p depends on z) gives us 


d^tu = --d^^p 
P 


dztw = -g 


'pd;,p-p'p\ pdzzp-p'd^p 


P^ 


Substituting these results into (IV.34) gives us 


P^ 


p' 


dttP = gpdtW + — dz:xP + gdzp - —p + dz:zP - —d^p 

P \ P P / 

1-0 ^ ( fi 9p' p' \ 

= —V P -f gpdtW + — gd^p - -p - —d^p 

P P \ P P ) 

Using (IV.31c) to remove the pdtw term from (IV.37) gives us 

o 1P^2 2 a ^ ^P ( a 9P' P' a \ 

OttP = —V p- g p- go^p + — gd^p - —p - —d^p 


P 


P 


P 


P 


IP 2 ( , 1PP'\ , IPf. ( . 1PP'\ ^ 

= yV U +^ I P + U +^ I 

If we solve (IV.31a) and (IV.31d) for dxU + dzW, we get 

gpw - dtp _ -p'w - dtp 
IP P 

which we can solve for w and differentiate with respect to time to get 

pdttP - ipdttp 


dtw = 


gp^ + 7pp' 


, P 7^ 0 


Equating this with dtW in (IV.31c) and solving for p gives 

P (pdttP - ipdttp 1 \ 

P = - - -2 , —, - -dzP 

g \ gp + ipp p J 


(IV.35) 

(IV.36) 

(IV.37) 

(IV.38) 

(IV.39) 

(IV.40) 

(IV.41) 
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Substituting this result into (IV.38) gives us 


dttP = 


f W -p{a + ^-f) 

[ +9fd,p-{g + af)d,p 

giPpdttp - gp^dttP IPp' ipdttp - pdttP 


= — V p — 
p 


gp"^ + 7PP' 


p 7P^ + 7pp' 


+ p^a,p(IV.42) 


after combining and cancelling terms. Now let us use an exponential atmospheric 
model. Combining (IV.32) and (11.59) gives us 


P = 


P = 


pa 

g 


pa 


Substituting these values into (IV.42) gives us 


dttP = —V p- 
P 

IP 


7P^2„ ai^^ttp - g^dttp 7y (idttp - ^dup" 


g^-iT 


g 


+ 


IP. 


\ 3$-Y- 


V p - '^-dttp + dttP + g—dzP 

pa p 

Noting that plp = g/a, this simplihes to 


+ g—o^p 
P 


(IV.43) 


dttp = V^p + gd^p 


(IV.44) 


We almost have a wave-like equation. The only thing missing is having the time 
derivative and the Laplace operator in the same variable. Look again at (IV.41). If 
we solve it for duP and use (IV.32) to cancel and combine terms, we have 

dttp = —dttP + ——— {gp + d^p) (IV.45) 

IP 7 

Substituting this result into (IV.44) gives us 

dttP = ? (V^p 4- gd^p) + p (1 - 7 ) pdtw , (IV.46) 

P ^ ' 

where we use (IV.31c) to make the substitution 


pdtW = -gp - d^p , 
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and we keep the usual ^ as the wave speed term rather than simplify it. If we try to 
remove the p and w terms, we go around in circles, never quite replacing them all in 
terms ofp. Hence, this is the closest we can come to deriving a wave-like equation that 
is equivalent to (IV.31). Our wave equation for p is altered by the vertical derivative 
of p and the time derivative of w. 

Note, however, that if we set 7 = 1, differentiate (IV.46) with respect to time, 
and make the appropriate substitutions, we can get 

dt (dttp - - PdztP = 0 , (IV.47) 

which more nearly resembles a wave equation, only for the time derivative of the 
pressure, rather than for the pressure itself. However, due to the dehnitions of the 
atmospheric constants, 7 = 1 means R = 0, which causes problems with the ideal 
gas law p = pRT on which we base our model. Hence, this contrivance is purely 
academic. 

3. Numerical Examples 

For our hrst example, we consider an open boundary on the sides of the do¬ 
main, so that the normal derivative is perpendicular to the force of gravity (see 
Fig. 19). We discretize (IV.31) according to the same second-order centered-difference 
scheme in Sec. B, which gives us 


\p = 

-p{A^u 

<1 

+ 

— p'w 

Atu = 

A^p 




P 



Atw = 

A^p 

9 

-P 



P 

P 


Atp = 

-IP ( 

<1 

+ 

s 

1 + gpw 


where we use z, w, p, and p instead of y, v, po, and po, respectively, to emphasize 
that our concern is now with a vertical, rather than horizontal, plane. As before, we 
use the same domain size and initial conditions. However, we now have hard walls 
on the top and bottom and NRBCs on the left and right, and the reference solution 
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Figure 19. An infinite channel domain truncated by artifical boundaries and 
Tr [After [22], Fig. 1, p. 2] 

is a 30 X 10 km channel, with the domain of interest in the center. The error norms 
are given in Table VI, with a comparison of pressure plots given by Fig. 20. 

Having successfully implemented gravity perpendicular to the open bound¬ 
aries’ normal vector, we now implement it parallel to the open boundary’s normal 
vector. Using again the same initial conditions, our domain is now a “bucket,” i.e., 
a semi-inhnite channel with the open boundary on top (Fig. 6). Table VII shows the 


J 

E, 

Eu 

Ejj 

Ep 

1 

0.23241 

0.52603 

0.12775 

0.23171 

2 

0.065394 

0.14839 

0.033298 

0.065251 

3 

0.03214 

0.069192 

0.017724 

0.032071 

4 

0.021193 

0.04438 

0.012525 

0.021144 

5 

0.015814 

0.031941 

0.0096649 

0.015776 

6 

0.012458 

0.024676 

0.007901 

0.012429 

7 

0.010184 

0.019914 

0.0066308 

0.01016 

8 

0.0086452 

0.016684 

0.0056677 

0.0086248 

9 

0.0076161 

0.01435 

0.0049472 

0.0075982 

10 

0.0068845 

0.012637 

0.0044394 

0.0068682 


Table VI. Error norms with J G 1... 10 in a horizontal channel with gravitational 
forces (IV.31) 
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Pressure perturbation. J=1 
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Pressure perturbation, J=10 



20 40 60 80 100 


Pressure perturbation error. J=1 



I 


Pressure perturbation error. J=10 


30 

20 

10 

0 

-10 

■20 

-30 



I 

I 


0,5 

0 

-0.5 

1 


Figure 20. Comparison of p in gravity-influenced system (IV.31) computed with J = 1 
and J = 10 in an infinite horizontal channel, with error norms computed by (IV. 14). 
(TL) Computed solution for J = 1. (TR) Computed solution for J = 10. (BL) 
Delta between reference solution and J = 1 computed solution. (BR) Delta between 
reference solution and J = 10 computed solution. 


error norms, and Fig. 21 compares the ^-velocity perturbations for J = 1 and J = 10. 

Next we combine these results into a domain open on three sides, with a hard 
wall on the bottom (see Fig. 22). Table VIII shows the error norms, and Fig. 23 
compares the x-velocity perturbations for J = 1 and J = 10. 
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J 

Ep 

Eu 

Ey 

Ep 

1 

0.097404 

0.10278 

0.22495 

0.096936 

2 

0.034899 

0.027336 

0.063934 

0.034738 

3 

0.026286 

0.014095 

0.030089 

0.026195 

4 

0.024522 

0.0098734 

0.019451 

0.024443 

5 

0.023927 

0.00764 

0.013899 

0.023852 

6 

0.023651 

0.0062309 

0.010768 

0.023578 

7 

0.023506 

0.0052081 

0.0086306 

0.023433 

8 

0.023425 

0.0044326 

0.0071758 

0.023353 

9 

0.023378 

0.0038462 

0.0060809 

0.023307 

10 

0.023348 

0.00339 

0.0052913 

0.023276 


Table VII. Error norms with J G 1... 10 in a vertical bncket with gravitational forces 
(IV.31) 


J 

E, 

Eu 

Eu 

Ep 

1 

0.31226 

0.85386 

0.30744 

0.30865 

2 

0.088812 

0.24972 

0.089607 

0.087862 

3 

0.043421 

0.11951 

0.04392 

0.042972 

4 

0.028607 

0.077987 

0.029393 

0.028308 

5 

0.021264 

0.057307 

0.021724 

0.021042 

6 

0.016764 

0.04499 

0.017264 

0.016589 

7 

0.01378 

0.036859 

0.014216 

0.013637 

8 

0.011743 

0.031518 

0.012194 

0.01162 

9 

0.010292 

0.027596 

0.010665 

0.010185 

10 

0.0091922 

0.024287 

0.0094066 

0.0090963 


Table VIII. Error norms with J G 1... 10 in an open-air domain with gravitational 
forces (IV.31) 
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2-velocity perturbation, J*1 
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Figure 21. Comparison of w in gravity-influenced system (IV.31) computed with 
J = 1 and J = 10 in a semi-infinite vertical channel, with error norms computed by 
(IV.14). (TL) Computed solution for J = 1. (TR) Computed solution for J = 10. 
(BL) Delta between reference solution and J = 1 computed solution. (BR) Delta 
between reference solution and J = 10 computed solution. 
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Figure 22. A half-plane domain D truncated by artihcial boundaries Fr, F^ and Fr 
[A fter [20], Fig. 1, p. 3] 
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Figure 23. Comparison of u in gravity-influenced system (IV.31) computed with 
J = 1 and J = 10 in an open-air domain, with error norms computed by (IV. 14). 
(TL) Computed solution for J = 1. (TR) Computed solution for J = 10. (BL) 
Delta between reference solution and J = 1 computed solution. (BR) Delta between 
reference solution and J = 10 computed solution. 
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E. ADVECTION 


Now we consider the case of non-zero mean flow. Beginning with the simple 
case, we have 


dtp + uod^p -t- vodyp + po {d^u + dyv) = 0 

dtu -f uodxU + vodyU H- d^p = 0 

Po 

dtv + u^dxV + vgdyV -I- dyp = 0 

Po 

dtp + uod^p + vodyP + 7Po {d^u + dyv) = 0 


(IV.49) 


Dehne a Lagrangian derivative 


Ttt — dt-\- uodx + vody 
Dt 


If we nse this dehnition in (IV.49), we get 

Dp 


— + Po {dxU + dyv) = 0 

Du 1 „ 

■!- ^xP = 0 

Dt Po 

Dv 1 ^ 

H- dyp = 0 

Dt Po 

^ + 7Po {dxU + dyv) = 0 , 


(IV.50) 


(IV.51) 


and we can easily show that this formnlation is eqnivalent to a Lagrangian scalar 
wave eqnation 


D^p _ 1P0^2 
DD po ^ 

If we incorporate Coriolis forces, then onr eqnation set is 

^ + Po {dxU dyv) = 0 

Du 1 . „ . . 

+ —dxP = f{v + Vo) 
Dt Po 

— + —dyP = -f{u + Uo) 
JJt Po 


^ + Wo {dxU + dyv) = 0 


(IV.52) 


(IV.53) 
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Using our now-familiar conversion technique, we have 

^ - 7Po/ {d^v - dyu) (IV.54) 

JJt Pq 

Now (IV.49) and (IV.53) both imply ^ so if we apply ^ to (IV.54), we 

get 

^ ~ ^ ~ ^ ° (IV.55) 

Using the same approach we used to convert the no-mean-flow form to the Klein- 
Gordon equation, we apply dy to (IV.53b), to (IV.53c), and subtract to get 

^ {da^v - dyu) = -f {da^u + dyv) (IV.56) 


Substituting this result into (IV.55) removes the vorticity term and gives us 

_D /DV ™ \ Cp _ 

Dt \Dt^^ Po ) Dt 

Integrate along characteristics, assuming an initial value of zero, to get 

D^p 

V = 

Po 


Dt^ 


p = — f^p . 


(IV.57) 


(IV.58) 


So our system is still equivalent to the Klein-Gordon equation, using Lagrangian 
derivatives rather than time derivative. However, we will show in Sec. 3 that such a 
system is inherently unstable under long-time integration. 


1. Interior Discretization Scheme 

Similar to (IV.50), if we dehne the Lagrangian difference 

A = At + uqAx + v^Ay , 

then the discrete basic system is given by 

Ap 4- Po {A^u -f Ayv) = 0 

V 1 . 

Am H-Aa;P = 0 

Po 

V 1 . 

Av H- Ayp = 0 

Po 

Ap + 7Po (Aa-M -4 Ayv) = 0 , 


(IV.59) 


(IV.60) 
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which is easily shown to be equivalent to 

AAp = — {A^A^p + AyAyp) (IV.61) 

Po 

Hence, the basic system under advection should be stable when implemented in con¬ 
junction with the Higdon NRBC dehned for a double-sized grid. Likewise, if we 
include Coriolis, we have the discrete system 


Ap + Po {A^u + Ayv) 

Au H- AxP 

Po 

Av + —AyP 

Po 

Ap + 7Po {^xU + Ayv) 


0 

/ + ^ + 

-/ (—p + U + Uo 
\Po 

0 , 


(IV.62) 


and if we use the discrete analog to the continuous case derivation in the preceding 
section, we can easily show that this system is equivalent to 

IPo 


AAp = — {A^A^P + Ay + Ayp) - fp 
Po 


(IV.63) 


which again should be stable when the Higdon scheme is used on a double-sized 
grid. Since the inclusion of gravity creates a system which is not compatible in the 
continuous case with a wave-like equation, its discrete formulation will also not be 
equivalent to a discrete wave-like equation. The discretization scheme in the xz plane 
under the influence of gravity is given by 

Ap + p [A^u + Azw) = -p' {w -f Wo) 

Au -|- -AxP = 0 
P 

V 1 A 9 

Aw +-AzP = --p 
P P 

Ap + 7p {A^u -f A^w) = gp{w^ wo) , (IV.64) 


where we replace the constants po and po with the ^-dependent values p and p. Note 
that Wo appears on the right-hand side of the equations for p and p because of the 
linearization process in Sec. H.C.2. In practice, however, the presence of the ground 
will require Wq = 0. 
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2. NRBC Formulation 

The addition of advection changes onr wave speed with respect to the bonnd- 


ary. To determine the new “correct” wave speed estimate, replace the time derivative 
in the Higdon NRBC dehnition (III. 11) with the Lagrangian derivative (expressed for 
a right-side bonndary): 


j 

n 


-1- Cjdx ) M = 0 

Dt ^ 


(IV.65) 


For waves striking the bonndary, we are only concerned with the velocity normal to the 
bonndary. Thns, we remove the tangential component of the Lagrangian derivative, 
expand and combine terms, and nse onr simplihcation Cj = cq, giving 


[dt + (co + Mo) dxY u = 0 


(IV.66) 


Since onr interior discretization scheme is eqnivalent (or approximately so) to a wave 
eqnation on a donble-sized grid, onr NRBC discretization for these eqnations is given 
by 


I - 
2St 


+ (co + Mo) 


2dx 


j 

a = 0 


(IV.67) 


3. Long-Term Instability of Advection with Coriolis 

Let ns briefly consider the discrete form of (11.63) in the xy plane with Coriolis 
forces. Using onr nsnal second-order centered-difference scheme, we have 


Atp + uqA^p + voAyp -F poA^u + poAyV 

Atu + uqA^u + voAyU + —A^p 

Po 

Atv -1- uqA^v + v^AyV H- Ayp 

Po 

Atp -f uqA^p + v^Ayp -f jpqA^u + jpoAyV 


0 

f{v + Vo) 

-/(mH-Mo) 

0 (IV.68) 


Assnme everything is initially at rest (no pertnrbations), and moMq ^ 0. With every 
pertnrbation variable eqnal to zero at times n < 0, onr valnes at time n = 1 become 



(IV.69) 


Ki 

2St 


= fvo 


m = 

P] 

^ = 0 
2St 

So, having started from a zero-perturbation domain, the combination of advection 
and Coriolis creates non-zero perturbations in u and v. In fact, u and v are still 
constant in the domain, but now the constant is non-zero. We note that p and p are 
still uniformly zero throughout the domains, and it is easy to show that they will 
always be zero. However, if we look at u and v at the next time step, noting that all 
spatial differences are zero, we get 


= -/K + «o 


(IV. 70) 


Substituting the previously determined values for u}, and nl , we get 


ulj = 2Stfvo-{2Stffuo 

vl = -2Stfuo-{2Stffvo (IV. 71) 

Proceeding to time step n = 3, we have 

ulj ^ 2 {2Stfvo) + O (^{2Stff^ 

~ -2 (25t/Mo) + O ((25t/)^) 

(IV. 72) 


where we are only concerned with the terms which are linear in 2Stf. It can be shown 
that 


]j ~ ( 25 ^/) vo + 0 [{2Stff^ 

j ^ -[^\{2Stf)uo + 0{{2Stff) , (IV.73) 
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where [-J denotes the floor function. 

Remark 1: If / 7 ^ 0 and either mq 7 ^ 0 or uq 7 ^ 0, then our zero-perturbation 
domain grows linearly, without bound, in u and/or v. While the growth will be slow 
due to the magnitude of /, it will still be present. Hence, the combination of Coriolis 
forces and non-zero advection is inherently unstable in the linearized Euler equations. 
We can only use this equation set for short-time integrations, not for longer durations. 

Remark 2: If we perform a similar analysis in the case of gravity and non¬ 
zero advection, we see that that system is inherently unstable if Wq 7 ^ 0. Horizontal 
advection does not cause any problems. Any physical problem which considers gravity 
and vertical advection is unsuited to this equation set. Meteorological problems, 
however, will generally include the ground and thus not have any constant vertical 
advection. 

4. Numerical Examples 

a. Basic System, Infinite Channel 

For the examples in this section, we use the same domain and initial 
conditions as before. Our first example is a basic system in an in fini te channel with a 
constant lateral advection of Uq = 100^. Again running the simulation up until t = 
24, we get the error norms listed in Table IX. Fig. 24 shows the density preturbation 
p for the J = 10 case. Note the faint wave propagating in the opposite direction; this 
false wave is a consequence of using centered-differences for our spatial discretization 
[24]. Fig. 25 compares the pressure perturbation p for low- and high-order NRBCs. 
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J 

E, 

Eu 


Ep 

1 

0.22447 

0.30496 

0.18074 

0.22447 

2 

0.06216 

0.09791 

0.059247 

0.062159 

3 

0.026666 

0.041939 

0.021853 

0.026666 

4 

0.017115 

0.025099 

0.012509 

0.017115 

5 

0.01258 

0.019265 

0.0097557 

0.01258 

6 

0.010007 

0.014986 

0.0074502 

0.010007 

7 

0.0082563 

0.012216 

0.0059657 

0.0082562 

8 

0.0070098 

0.010347 

0.0051175 

0.0070098 

9 

0.0060805 

0.0090438 

0.0044283 

0.0060805 

10 

0.0053682 

0.0079242 

0.003876 

0.0053683 


Table IX. Error norms with J G 1... 10 in an infinite channel with horizontal advec- 
tion 



Fignre 24. Plot of p in basic system (IV.49) with left-to-right advection with J = 10 
in an infinite channel. (BL) Compnted solntion. (Top) Reference solntion; the area 
corresponding to the compnted solntion is contained between the vertical lines. (BC) 
Reference solntion tru n cated to compnted solntion domain. (BR) Delta between 
reference solntion and compnted solntion, with error norm compnted by (IV. 14). 
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Figure 25. Comparison of p in basic system (IV.49) with left-to-right advection 
computed with J = 1 and J = 10 in an in fini te channel, with error norms computed 
by (IV.14). (TL) Computed solution for J = 1. (TR) Computed solution for J = 10. 
(BL) Delta between reference solution and J = 1 computed solution. (BR) Delta 
between reference solution and J = 10 computed solution. 
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J 

E, 

Eu 

Ey 

Ep 

1 

0.33816 

0.36566 

0.36566 

0.33819 

2 

0.091059 

0.10563 

0.10563 

0.091066 

3 

0.043514 

0.050827 

0.050827 

0.043517 

4 

0.028633 

0.032601 

0.032601 

0.028635 

5 

0.021446 

0.023794 

0.023794 

0.021447 

6 

0.0172 

0.018581 

0.018581 

0.017201 

7 

0.014495 

0.015508 

0.015508 

0.014496 

8 

0.012603 

0.013347 

0.013347 

0.012604 

9 

0.011967 

0.013011 

0.012734 

0.011969 

10 

0.078775 

0.10216 

0.082933 

0.078757 


Table X. Error norms with J G 1... 10 in an open domain with diagonal advection 


b. Basic System, Open Domain 

Now consider the basic system in an open domain with a constant 
diagonal advection of Uq = 90 m/s, vq = 90 m/s. This time, we get the error norms 
listed in Table X. Not snrprisingly, the error norms for u and v are identical np to 
J = 8 (with V slightly better than u in the J = 9 case). In this example, we see that 
the J = 10 error norms are larger than those for J = 9. Fig. 26 shows the a:-velocity 
pretnrbation u for the J = 8 case; Figs. 27 and 28 shows the same variable for the 
J = 9 and J = 10 cases, respectively. We see from the latter two that the system 
begins to destabilize at the top-right corner. The reason for this destabilization is 
still nnder investigation. Fig. 29 compares the |/-velocity pertnrbation v for low- and 
high-order NRBCs. 
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Error nomv-0 010936 
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Figure 26. Plot of u in basic system (IV.49) with bottom-left-to-top-right advection 
with J = 8 in an open domain. (TL) Computed solution. (Right) Reference solution; 
the area corresponding to the computed solution is contained in the center box. (CL) 
Reference solution truncated to computed solution domain. (BL) Delta between 
reference solution and computed solution, with error norm computed by (IV. 14). 
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Error t>omv-0.009933$ 


Figure 27. Plot of u in basic system (IV.49) with bottom-left-to-top-right advection 
with J = 9 in an open domain. (TL) Computed solution. (Right) Reference solution; 
the area corresponding to the computed solution is contained in the center box. (CL) 
Reference solution truncated to computed solution domain. (BL) Delta between 
reference solution and computed solution, with error norm computed by (IV. 14). 
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Figure 28. Plot of u in basic system (IV.49) with bottom-left-to-top-right advection 
with J = 10 in an open domain. (TL) Computed solution. (Right) Reference solution; 
the area corresponding to the computed solution is contained in the center box. (CL) 
Reference solution truncated to computed solution domain. (BL) Delta between 
reference solution and computed solution, with error norm computed by (IV. 14). 


Y-ve*odty perturtaUon, J*1 Y-velocity perturbation. J=6 
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Figure 29. Comparison of v in basic system (IV.49) with bottom-left-to-top-right 
advection computed with J = 1 and J = 8 in an open domain, with error norms 
computed by (IV. 14). (TL) Computed solution for J = 1. (TR) Computed solution 
for J = 8. (BL) Delta between reference solution and J = 1 computed solution. (BR) 
Delta between reference solution and J = 8 computed solution. 
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J 

E, 

Eu 

Ey 

Ep 

1 

0.22154 

0.59249 

0.25232 

0.22154 

2 

0.10104 

0.3566 

0.15752 

0.10104 

3 

0.34963 

1.3304 

0.39134 

0.34963 

4 

1.4844 

5.1058 

1.391 

1.4844 

5 

2.3286 

7.4879 

3.3526 

2.3286 

6 

34.219 

98.952 

36.169 

34.219 

7 

71.57 

277.37 

83.65 

71.57 

8 

4971.8 

14265 

4154 

4971.8 


Table XL Error Norms with J G 1... 8 in an In fini te Channel with Horizontal Ad- 
vection and Coriolis 


c. Coriolis Forces, Infinite Channel 

Having shown already that the system is inherently nnstable for long¬ 
time integrations involving both advection and Coriolis, let ns look briefly at the 
effectiveness of the Higdon NRBCs for a short-dnration simnlation. We first consider 
an in fini te channel with left-to-right advection, Uq = 100^. Incorporating Coriolis 
forces, we nse the discretization scheme (IV.62) for the interior. Table XI shows the 
error norms for J G 1... 8. We see some improvement between J = 1 and J = 2, bnt 
the system fails for higher J. We will not show any plots of these resnlts. 
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J 

E, 

Eu 

Ey 

Ep 

1 

0.86499 

0.31392 

0.33843 

0.86506 

2 

0.58828 

0.26882 

0.2833 

0.58833 

3 

6.288 

2.5314 

2.67 

6.2885 

4 

327.16 

134.55 

141.92 

327.18 

5 

14895 

4854.5 

5120.1 

14896 

6 

278880 

87044 

91822 

278900 


Table XII. Error Norms with J G 1... 6 in an Open Domain with Diagonal Advection 
and Coriolis 


d. Coriolis Forces, Open Domain 

Likewise, when we consider the system in an open domain, with the 
NRBC on all fonr sides, nsing a diagonal advection, uq = Vq = 100^, we get the 
resnlts shown in Table XII. A slight improvement from J = 1 to J = 2 is followed by 
massive instability for higher J. Again, we will omit any plots of these resnlts. 
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J 

E, 

Eu 


Ep 

1 

0.2251 

0.29999 

0.18206 

0.22442 

2 

0.062733 

0.096051 

0.059583 

0.062517 

3 

0.027088 

0.041293 

0.022049 

0.027048 

4 

0.017524 

0.024613 

0.012558 

0.017488 

5 

0.012913 

0.018886 

0.0098005 

0.012887 

6 

0.010289 

0.014672 

0.0074523 

0.01027 

7 

0.0085063 

0.011944 

0.0059542 

0.0084902 

8 

0.007226 

0.010119 

0.005107 

0.0072123 

9 

0.0062668 

0.0088368 

0.0044166 

0.0062548 

10 

0.0055272 

0.0077406 

0.0038622 

0.0055167 


Table XIII. Error norms with J G 1... 10 in an infinite channel with horizontal 
advection and gravity 


e. Gravity, Infinite Channel 

Tnrning onr attention now to the inclnsion of gravity, we again consider 
an in fini te channel with horizontal advection, uq = 100^. Using the discretization 
scheme (IV.64), we get the error norms shown in Table XIII. Unlike the Coriolis 
case, this time we get good resnlts all the way np to J = 10. Fig. 30 illnstrates the 
difference in pressnre pertnrbation resnlts for small and large J. 
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Pressure perturbation. J=1 
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Figure 30. Comparison of p in gravity-influenced system (IV.64) with left-to-right 
advection computed with J = 1 and J = 10 in an inhnite channel, with error norms 
computed by (IV. 14). (TL) Computed solution for J = 1. (TR) Computed solution 
for J = 10. (BL) Delta between reference solution and J = 1 computed solution. 
(BR) Delta between reference solution and J = 10 computed solution. 
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J 

E, 

Eu 


Ep 

1 

0.25917 

0.35427 

0.28085 

0.25639 

2 

0.072395 

0.11036 

0.084279 

0.071655 

3 

0.032031 

0.050496 

0.036455 

0.03179 

4 

0.020847 

0.03178 

0.022262 

0.020674 

5 

0.015383 

0.023689 

0.016319 

0.015259 

6 

0.012198 

0.018435 

0.012387 

0.012104 

7 

0.010073 

0.01515 

0.0099833 

0.0099939 

8 

0.0085547 

0.012836 

0.0084265 

0.0084879 

9 

0.0074187 

0.011229 

0.007244 

0.0073608 

10 

0.0066163 

0.010216 

0.0066934 

0.0065656 


Table XIV. Error norms with J G 1... 10 in a half-plane with horizontal advection 
and gravity 


/. Gravity, Ground and Open Air 

Looking now at an “open” domain in the presence of gravity, we con¬ 
sider the open-air domain of Fig. 22. This domain matches the physical idea of 
modeling an open atmosphere bonnded by the gronnd. While the presence of the 
gronnd prevents ns from considering vertical or diagonal advection, it still presents 
ns with a domain containing two adjacent open bonndaries. We keep uq = 100^ for 
onr advection speed and rnn onr standard simnlation. The error norms for different 
valnes of J are given in Table XIV. Fig. 31 illnstrates the differences between low 
and high J. 

F. SUMMARY 

This high-order NRBC implementation provides high accnracy with little com- 
pntational overhead. However, it has three signihcant limitations: 

1. The formnlation reqnires one-sided high-order spatial and temporal deriva¬ 
tives. 

2. These one-sided derivatives limit the NRBC order to the size of the domain. 

3. Coriolis and advection cannot be nsed together stably. 
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Figure 31. Comparison of p in gravity-influenced system (IV.64) with left-to-right 
advection computed with J = 1 and J = 10 in an open-air domain, with error norms 
computed by (IV. 14). (TL) Computed solution for J = 1. (TR) Computed solution 
for J = 10. (BL) Delta between reference solution and J = 1 computed solution. 
(BR) Delta between reference solution and J = 10 computed solution. 


In the next chapter, we will look at another NRBC implementation, one which does 
not suffer from these limitations. 
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V. GIVOLI-NETA NRBCS FOR THE 
LINEARIZED EULER EQUATIONS 


A. INITIAL IMPLEMENTATION FOR FIRST-ORDER 
SYSTEMS 
1. Derivation 

Let us now consider the Givoli-Neta auxiliary variable NRBC described in 
Sec. III.B.2. That section described this NRBC formulation for a scalar wave-like 
equation. In this chapter we apply it to the linearized Euler system. We hrst dehne 
the auxiliary variables in vector form: 

(pj = (pj+i (V.l) 

(po = (p 


Note that we have left off the truncation condition (pj = 0 from our formulation. 
As we derive the implementation below, we shall see that we need to modify the 
truncation condition slightly; the modihed truncation condition will be given as part 
of the characteristic-based derivation (V.IO). 

Givoli and Neta [39, 40] show how to manipulate the Klein-Gordon equation to 
remove the normal derivatives. We will show briefly how to do so with the linearized 
Euler equations. For dehniteness, we proceed for the right-side boundary; thus, n = x. 

Begin with the equation system in vector form: 

dt0 + Adx0 + Bdyip = C(p + D , (V.2) 


where 


A = ( 

/ 


\ 

T 

kP 

U V 

p j 



' Mo 

Po 

0 

0 ^ 


0 

Uq 

0 

J_ 

A = 



Po 


0 

0 

Mo 

0 


1 0 

IPo 

0 

Mo ) 
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^ 1^0 0 Po 0 ^ 

0 Vq 0 0 

B = 

0 0 Uo - 

^ po 

0 0 7Po ^'o / 

In the basic system, C and D are both zero. With Coriolis, we have 

^ 0 0 0 0 ^ 

0 0/0 
0-/00 
^ 0 0 0 0 y 

D = 0 /no -fuo 0 ) 

With gravity, replace v, Vq, Pq, and po in A and B with w, Wq, p, and p, respectively, 
and use 


C = 


D = 


0 0 -p' 0 

0 0 0 0 

-p4 0 0 0 

0 0 pp 0 y 

-p'wo 0 0 gpwo 


Dehne the linear differential operator L{(p) to be 


L{0) = dtp + AdxP + Bdyp — C(p (V.3) 

By dehnition, I/((p) = D. A more useful result is the following lemma: 

Lemma V.l The auxiliary variables pj satisfy 

i(ft) = 0 (V.4) 

for all j G 1... J. The only exception to this lemma is when both of the following 
conditions are met: 
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• n = i/c 

• ^ 0 

Proof. Apply the operator L to (V.l) to get 


1 


L{ipi) = L dH + —dt (p 


1 


= dH{L{<p)) +—dt{L{^)) 
= dfjD + -d,D 


= 0 


(V.5) 


Then proceed by induction to prove the lemma for all j □. 

(Note that if both conditions in Lemma V.l’s exception are met, then L {dz<p) ^ 
dz {L and we cannot proceed with this proof or the overall derivation. The NRBC 
which handles these conditions will be given in Sec. C.) 

Solve (V.l) for dz;<Pj and, using Lemma V.l, substitute the result into (V.2): 


Aipj+i = —A - I dt<Pj - Bdt(pj + CLpj , 'ij el... J 


(V.6) 


We now have an auxiliary variable equation which need be dehned only on the bound¬ 
ary. However, as it is dehned solely on the boundary, and given that everything on 
the boundary is initially zero, we have an equation system which will always be iden¬ 
tically zero. We will take a characteristic-based approach similar to that presented 
in Sec. 6 of [57]. 

Let A be the diagonal matrix of the eigenvalues of A in decreasing order, with R 
the corresponding right eigenvectors. With two degrees of freedom for the eigenvectors 
associated with A = mq; we choose eigenvectors which do not map one characteristic 

variable to one state variable, which we would get if we chose ^ 1 0 0 0 ^ and 

W 

0 0 10 for our uq eigenvectors. Instead, we use the following: 


A = diag ^ tiQ + Co Uq Uq Uq- 
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R 

R-^ 

A 


PO I I 

CO 

1 0 0 

0 -1 1 

poco 0 0 




(V.7) 


With this dehnition, left-multiply (V.6) by R ^ and insert RR ^ as needed to get 


A6 

A^+i 



d/o - m{o + cco + n 

dtCj - mCj + CCj, , 


(V.8) 


where 


0 Aj 

B = R-^BR 
C = R-^CR 
D = R~^D 


(V.9) 


However, we still have a system which is dehned exclnsively on the bonndary. To over¬ 
come this problem, nse (V.6) in characteristic form as a natnral bonndary condition 
for the ontgoing characteristics. To close the system, nse a Dirichlet condition for the 
highest-order anxiliary variables corresponding to the incoming characteristics. The 
hnal system then is 

dt^o = —Adx^o — Bdy^o + + D ^ 

= —Bdy^o -f- C^o + D 
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/--A = -BdyC, + CCj , j el...J 


0j = 0 


(V.IO) 


where the subscripts + and — mean we apply these equations only to the outgoing 
and incoming characteristics, respectively; that is, the hrst equation only applies to 
the outgoing characteristic variables, the second and third equations apply to all four 
characteristic auxiliary variables (for each j), and the hnal equation only applies to 
the incoming characteristic auxiliary variables. For example, if we have J = 1, Mq > 0, 
C = 0, D = 0, and we use standard hrst-order forward-differences in time, hrst-order 
backward-differences in x, and second-order centered-differences in y, then (V.IO) 
dehnes the matrix system 


1 

0 

0 

0 

0 

0 

0 

0 


dn 

^a,E,i 

0 

1 

0 

0 

0 

0 

0 

0 


tn 

^b,E,i 

0 

0 

1 

0 

0 

0 

0 

0 


tn 

^c,E,i 

1-A„ 

0 

0 

0 

StXa 

0 

0 

0 


^d,E,i 

0 

1 — Af, 

0 

0 

0 

StXb 

0 

0 


tn 

^l,a,E,i 

0 

0 

1-Ae 

0 

0 

0 

StXc 

0 



0 

0 

0 

1 — Arf 

0 

0 

0 

StXd 


dn 

^l,c,E,i 

. 0 

0 

0 

0 

0 

0 

0 

1 ; 


tn 

\ ^l,d,EA / 


/ en-i 
^a,E,i 




n—1 

b,E,i 

cn—1 

'!>c,E,i 



(V.ll) 
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where A denotes the individual eigenvalues of the diagonal matrix A; the subscript 
1 denotes the hrst (and only) auxiliary variable; the subscripts a, b, c, d denote the 
individual components of the characteristic/state/auxiliary variable vectors; the sub¬ 
scripts E and E — 1 denote the eastern boundary and the point adjacent to it, respec¬ 
tively; the subscripts i, i -|- 1, and i — 1 denote the grid point in the y direction; and 
the superscripts n and n — 1 denote the current and previous time steps, respectively. 

Note that if we are applying this NRBC in the presence of gravity, even in a 
horizontal channel, we must take care to consider the ^-derivative of the eigenvector 
matrix R when we convert from state variables to characteristics. In such a case, 
(V.IO) becomes 




Ci / 



dtij -l- A^_|_i 

^3 / 


-Ad^fo - + {b R + C) + D 

-Bd^io + {b (d,R-^) R + C)i, + D 

-ml A {b (d,R-^) R + c)i, , iel...J 


0J = o] _ . (V.12) 

2. Incompatibility with Zero Advection 

Note that, because of the A on the left-hand side of (V.10b,c), all of the 
eigenvalues must be non-zero. Requiring the advection be subsonic is trivial; it can 
be a constraint of the model. But the case of zero advection is trickier. The answer is 
to dehne a “false” non-zero advection, small enough that there will be no drift from 
one point to an adjacent point over the duration of the simulation: 

min(5a:) 


Mo < 


2A 


(V.13) 


where min(5a:) is the smallest grid spacing in the x direction, and tmax is the duration 
of the simulation. If the false Uq is too small, however, there is a risk that the resulting 
matrix could be numerically singular. Realistically, this concern is unfounded; we have 
seen experimentally (using Matlab) that the matrix A is not ill-conditioned so long 
as ||mo|| > (approximately an eighth of an inch per year). 



J 

E, 

Eu 

Ey 

Ep 

1 

0.15985 

0.32103 

0.13114 

0.15985 

2 

0.1581 

0.32249 

0.12376 

0.1581 

3 

0.15697 

0.32307 

0.12029 

0.15697 

4 

0.15617 

0.32308 

0.11844 

0.15617 

5 

0.15557 

0.32274 

0.11729 

0.15557 

6 

0.15506 

0.32218 

0.11646 

0.15506 

7 

0.15462 

0.32151 

0.11579 

0.15462 

8 

0.15422 

0.32077 

0.1152 

0.15422 

9 

0.15384 

0.32 

0.11465 

0.15384 

10 

0.15347 

0.3192 

0.11413 

0.15347 


Table XV. Error norms for Givoli-Neta NRBCs in basic system (IV.1) with J G 
1... 10 in an infinite channel with no advection 


3. Numerical Examples 

Let ns consider some nnmerical examples. In each example, we nse onr stan¬ 
dard domain and initial condition in an inf inite channel (Fig. 19). We hrst consider 
the basic case with no advection. Since we have to dehne a false advection on the 
bonndary, we set 

uo = (V.14) 

^^max 

For this domain and simnlation dnration, that eqnates to Uq = 1.3889^. Fnrther- 
more, we set Cj = 1 for all j [23, 55], a simplihcation we will nse for all snbseqnent 
derivations for these NRBCs. Table XV shows the state variable error norms for 
J G 1... 10, Fig. 32 shows the state variable p for the J = 10 case, and Fig. 33 com¬ 
pares the state variable u for the J = 1 and J = 10 cases. We note that there is very 
little difference between the low-order and high-order resnlts. This is a conseqnence of 
the characteristic-based approach [72] as well as a conseqnence of converting the nor¬ 
mal derivative to tangential [36]. If we instead nse a horizontal advection uq = 100^, 
we get the resnlts shown in Table XVI and Figs. 34 and 35. 

If we incorporate Coriolis forces in this domain, we see that the Givoli-Neta 
NRBCs still prodnce valid resnlts. Tables XVII and XVIII show the error norms for 
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Figure 32. Plot of p for Givoli-Neta NRBCs in basic system (IV. 1) with J = 10 in 
an infinite channel with no advection. (BL) Computed solution. (Top) Reference 
solution; the area corresponding to the computed solution is contained between the 
vertical lines. (BC) Reference solution truncated to computed solution domain. (BR) 
Delta between reference solution and computed solution, with error norm computed 
by (IV. 14). 


the state variables for the Mq = 0 and Uq = 100 cases, respectively. In the interest of 
brevity, we will omit any hgures from these examples. 

Similarly, incorporating gravity in the xz plane poses no additional difficulties 
in this horizontal channel, as shown in Tables XIX and XX. We will see in Sec. C 
how to deal with the effects of gravity on a vertical open boundary. 

Interestingly, even though these error norms appear to decrease only slightly 
as J increases, we can see from a logarithmic plot that they do, in fact, decay ex¬ 
ponentially. See Fig. 36 for a graphical representation of the values in Table XX. 
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J 

E, 

Eu 


Ep 

1 

0.1556 

0.23218 

0.15239 

0.1556 

2 

0.15506 

0.23139 

0.15096 

0.15506 

3 

0.15452 

0.2306 

0.14954 

0.15452 

4 

0.15399 

0.22982 

0.14815 

0.15399 

5 

0.15346 

0.22904 

0.14679 

0.15346 

6 

0.15293 

0.22827 

0.14545 

0.15293 

7 

0.1524 

0.22751 

0.14413 

0.1524 

8 

0.15188 

0.22674 

0.14283 

0.15188 

9 

0.15136 

0.22599 

0.14156 

0.15136 

10 

0.15085 

0.22523 

0.1403 

0.15085 


Table XVI. Error norms for Givoli-Neta NRBCs in basic system (IV.1) with J G 
1... 10 in an infinite channel with left-to-right advection 


J 

E, 

Eu 

Eu 

Ep 

1 

0.12541 

0.25444 

0.1244 

0.12541 

2 

0.12484 

0.25298 

0.11594 

0.12484 

3 

0.12438 

0.25176 

0.11182 

0.12438 

4 

0.12394 

0.25063 

0.10959 

0.12394 

5 

0.12351 

0.24953 

0.10822 

0.12351 

6 

0.12307 

0.24843 

0.10724 

0.12307 

7 

0.12263 

0.24734 

0.10647 

0.12263 

8 

0.12218 

0.24626 

0.10581 

0.12218 

9 

0.12174 

0.24518 

0.10521 

0.12174 

10 

0.1213 

0.24411 

0.10465 

0.1213 


Table XVII. Error norms for Givoli-Neta NRBGs in Goriolis-inflnenced system (IV. 15) 
with J G 1... 10 in an inf ini te channel with no advection 
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J 

E, 

Eu 

Ey 

Ep 

1 

0.24986 

0.6637 

0.26511 

0.24986 

2 

0.24943 

0.66284 

0.26447 

0.24943 

3 

0.24901 

0.662 

0.26385 

0.24901 

4 

0.2486 

0.66119 

0.26325 

0.2486 

5 

0.24819 

0.66039 

0.26268 

0.24819 

6 

0.24779 

0.65961 

0.26213 

0.24779 

7 

0.2474 

0.65886 

0.2616 

0.2474 

8 

0.24702 

0.65812 

0.26109 

0.24702 

9 

0.24664 

0.6574 

0.2606 

0.24664 

10 

0.24626 

0.65669 

0.26013 

0.24627 


Table XVIII. Error norms for Givoli-Neta NRBCs in Coriolis-inflnenced system 
(IV. 15) with J G 1... 10 in an infinite channel with left-to-right advection 


J 

E, 

Eu 

Eu 

Ep 

1 

0.12998 

0.25364 

0.1265 

0.1294 

2 

0.1294 

0.25218 

0.11796 

0.12881 

3 

0.12893 

0.25098 

0.11379 

0.12833 

4 

0.12849 

0.24986 

0.11152 

0.12788 

5 

0.12804 

0.24877 

0.11012 

0.12742 

6 

0.12758 

0.24769 

0.10913 

0.12697 

7 

0.12713 

0.24662 

0.10835 

0.12651 

8 

0.12667 

0.24555 

0.10768 

0.12605 

9 

0.12621 

0.24449 

0.10707 

0.1256 

10 

0.12575 

0.24343 

0.1065 

0.12514 


Table XIX. Error norms for Givoli-Neta NRBGs in gravity-inflnenced system (IV.31) 
with J G 1... 10 in an inf ini te channel with no advection 
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X-vek>city perturbation, J=1 



X-velocity perturbation, J=10 





Figure 33. Comparison of u for Givoli-Neta NRBCs in basic system (IV.1) computed 
with J = 1 and J = 10 in an infinite channel with no advection, with error norms 
computed by (IV. 14). (TL) Computed solution for J = 1. (TR) Computed solution 
for J = 10. (BL) Delta between reference solution and J = 1 computed solution. 
(BR) Delta between reference solution and J = 10 computed solution. 


J 

E, 

Eu 


Ep 

1 

0.1402 

0.20569 

0.14568 

0.13973 

2 

0.13959 

0.20481 

0.14417 

0.13912 

3 

0.13899 

0.20395 

0.14269 

0.13852 

4 

0.13839 

0.20309 

0.14123 

0.13793 

5 

0.13779 

0.20223 

0.13979 

0.13733 

6 

0.1372 

0.20138 

0.13837 

0.13674 

7 

0.1366 

0.20053 

0.13698 

0.13615 

8 

0.13601 

0.19969 

0.13561 

0.13556 

9 

0.13543 

0.19885 

0.13426 

0.13498 

10 

0.13484 

0.19802 

0.13293 

0.1344 


Table XX. Error norms for Givoli-Neta NRBCs in gravity-influenced system (IV.31) 
with J G 1... 10 in an inhnite channel with left-to-right advection 
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X-direction velocity perturtation. reference solution 



•50 0 50 100 150 200 

X-direction velocity perturbation Truncated reference solution Solution delta 




50 100 SO 100 50 100 

Error norm-0.19874 



Figure 34. Plot of u for Givoli-Neta NRBCs in basic system (IV.1) with J = 10 
in an infinite channel with left-to-right advection. (BL) Computed solution. (Top) 
Reference solution; the area corresponding to the computed solution is contained 
between the vertical lines. (BC) Reference solution truncated to computed solution 
domain. (BR) Delta between reference solution and computed solution, with error 
norm computed by (IV. 14). 
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Figure 35. Comparison of p for Givoli-Neta NRBCs in basic system (IV. 1) computed 
with J = 1 and J = 10 in an infinite channel with left-to-right advection, with error 
norms computed by (IV. 14). (TL) Computed solution for J = 1. (TR) Computed 
solution for J = 10. (BL) Delta between reference solution and J = 1 computed 
solution. (BR) Delta between reference solution and J = 10 computed solution. 


Density perturbation error norm 


X-velocity perturbation error norm 




Z-velocity perturbation error norm Pressure perturbation error norm 




Figure 36. Logarithmic plot of state variable error norms (IV. 14) for Givoli-Neta 
NRBCs applied to gravity-influenced system (IV.31) with J G 1... 10 in an inhnite 
channel with left-to-right advection. (TL) Error norms for p. (TR) Error norms for 
u. (BL) Error norms for w. (BR) Error norms for p. 
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B. CORNER CONDITIONS IN AN OPEN DOMAIN 
1. Derivation 

In an anxiliary variable NRBC scheme, the corner where two open bonndaries 
intersect is a sonrce of concern. Assnming the tangential derivatives are approximated 
with a centered-difference formula, points adjacent to the corner depend on the value 
at the corner. (This is not the case in a standard high-order NRBC such as the 
Higdon sequence, where there are no tangential derivatives and thus no dependence 
on the corner values, as we noted in Sec. IV.B.4.) So how should we treat this corner? 
Is it to be considered part of one boundary but not the other? Even so, we still have 
to use a different discretization scheme than that used elsewhere on the boundary. 
We propose here to treat the corner as if it were part of a curved boundary; dehne the 
normal derivative as the 45° vector bisecting the normal vector of the two adjacent 
sides (see Fig. 37). With this dehnition, we have for the top-left and top-right corners 



Figure 37. The “normal” derivative for the top-right corner of an open domain 


of Fig. 22, 


Top-left: dn 
Top-right: dn 


-dx + d^ 

V2 

dx + d^ 

V2 ’ 


(V.15) 
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and the associated tangential derivatives are given by 


Top-left: dr 
Top-right: dr 


dr + d^ 

72 
dx - 
72 


(V.16) 


In order to derive onr characteristic system for the anxiliary variables, we mnst rewrite 
(V.2) in terms of normal and tangential derivatives. It is straightforward to show that 


dt(p + Ndn0 + Tdr(p = C(p -f D , 


(V.17) 


where 


Top-left: N 
T 

Top-right: N 
T 


1 

72 

1 

72 

1 

72 

1 


i-A + B) 
(A + B) 
(A + B) 
(A-B) . 


(V.18) 


With these dehnitions, it is easy to show that the characteristic NRBC system for 
each corner is given by 


'dtCo = -Ad,Co-fdrCo + CCo + D]^ 

(/-A)at^+Ag+1 = -fdrCj + CC, (V.19) 

[6 = o]_, 

where 

A = diag ^ no -f Co no no no - Cq ^ 

A = RNR-^ 

f = R-^TR , (V.20) 


and C and D dehned as in (V.9). On the top-left corner, we have 



To = 


R-^ = 


Uo + Wo 

V2 

( 1 


R = 


1 -1 


CQ 


\/2po 

_CQ_ 


\/2po 


CQ 

\/2po 

CQ 

\/2po 


0 0 

\/2po 


V 


4co 

1 

4 

1 

4 

\/2po 

4co 


\/2po 

4co 

1 

4 

1 

4 

_ \/2po 
4co 


2c; 


2c; 


and on the top-right corner, we nse 


no 

t-q 






Wo + Wo 

72 

Wo - Wo 


72 


/ 


CQ 

\/2po 

CQ 

\/2po 


''0 


-1 


n \/2po 
4co 


1 

4 

_ 1 
4 

\/2po 

4co 


1 1 \ 

_1 CQ 

\/2po 

1 CQ 

\/2po 

0 eg ^ 

\/2po 1 \ 

4co 2c^ I 


1 1 



\/2po _1 I 

4co ^ / 


In implementing the hnite-difference approximations of 5^-, decompose it into its con- 
stitnent dx and components. Use the hr anxiliary variable for dx-, and nse the F^ 
or Fr anxiliary variable for d^-, nsing one-sided differences. Note also that the eigen- 
valnes of N mnst be non-zero. We can again nse a “false” advection to get aronnd 
this difficnlty, if needed. 
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J 

E, 

Eu 

Ey 

Ep 

1 

0.16986 

0.20783 

0.19382 

0.16986 

2 

0.16904 

0.20686 

0.19168 

0.16905 

3 

0.16823 

0.20591 

0.18957 

0.16823 

4 

0.16742 

0.20496 

0.18749 

0.16742 

5 

0.16661 

0.20401 

0.18545 

0.16661 

6 

0.16581 

0.20307 

0.18344 

0.16581 

7 

0.16501 

0.20214 

0.18146 

0.16502 

8 

0.16422 

0.20121 

0.17951 

0.16422 

9 

0.16343 

0.20029 

0.1776 

0.16343 

10 

0.16265 

0.19938 

0.17571 

0.16265 


Table XXL Error norms for Givoli-Neta NRBCs in basic system (IV.1) with J G 
1... 10 in an open half-plane with left-to-right advection 

2. Numerical Examples 

For simplicity, let ns consider jnst one example of this open-air domain. We 
nse onr standard example with left-to-right advection with no inhomogeneons forces. 
We reqnire a “false” vertical advection on the top bonndary and at the corners; we 
nse Wq = = 1.0417^. Table XXI and Fig. 38 show the error norms for the 

state variables for J G 1... 10. Again, althongh the differences are small, they are 
nonetheless exponential. 

C. GRAVITATIONAL EFFECTS 
1. Derivation 

Onr claim that L{(pj) = 0 in the preceding section hinges on the fact that the 
coefficient matrices are constant in the direction of the normal derivative. However, 
this is not the case when we consider an NRBC on Ft in the presence of gravity. In 
that case, the coefficient terms A, B, C, and D all depend on Recall the differential 
operator L{0) dehned in (V.3). We still have L{(p) = D hy dehnition, bnt what of 
the anxiliary variables? Since L {8^0) 0 {L {0)), we come np against the exception 
to Lemma V.l. A new derivation is needed. Apply the operator L to both sides of 
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Figure 38. Logarithmic plot of state variable error norms (IV. 14) for Givoli-Neta 
NRBCs applied to basic system (IV.1) with J G 1... 10 in an open half-plane with 
left-to-right advection. (TL) Error norms for p. (TR) Error norms for u. (BL) Error 
norms for v. (BR) Error norms for p. 


(V.la, i = 0), and simplify 


i(w) 


i(^i) 


L(d,(^ + L(d,0) 

/ 

\dttp + Ad^t<^ - Cdt(p\ 

< 

+ [dzS + Ad^^(p - Cd:,(p] 

/ 

dt (dtp + Ad^p + Bd^p - Cp) 

< +d^ (dtp -f Ad^p + Bd^p - Cp) 

-(d^A)d^p- (d^B)d^p + (d^C)p 
dt(L(p)) + d^(L(p)) - (d^A)d^p- (d^B)d^p + (d^C)p 

dtD + d^D - (d^A)d^p- (d^B)dzp + (d^C)p 

-(d^A)da;p- (d^B)d^p + (d^C)p + d^D (V.21) 


As we continue applying the operator L to successive Pj's, this expression will result 
in ever-higher 2 ;-derivative terms, exactly the case we are trying to avoid. We need 
to dehne our reference states in such a way as to remove this difficulty. The most 
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effective approach is to use the exponentially-decaying atmospheric model (IV.32), 
given in Sec. IV.D.l. 

We can rewrite the matrices A and B as the sum of a constant matrix and a 
^-dependent one: 


A — rto/ + Az 

B = wqI Bz , (V.22) 

where Az and Bz are the ^-dependent cells of A and B, respectively. Given these 
^-dependencies, and the dehnitions for p and p, we have 


dzA = —aAz 
dzB = —aBz 
dzC = -aC 
dzD = -aD . 


Continuing our derivation of L{(pi) using these dehnitions, we have 


L{ipi) = aAzdx0 + aBzdz(p — aC(p — aD 

{l{(p) - D - dt(p - UQdx(p - WQdz(p 
{p - D - dtif - uodx<f - wodz(p 
= -a {dt(p + UQdx(p + WQdz(p) 

L{(pi) = -adpif , 


= a 


= a 


(V.23) 


(V.24) 


where dp denotes the Lagrangian derivative along the how. If we continue this deriva¬ 
tion for successive (pj terms, we get 

L{02) = —2adp(pi — a^dpip 

L{(p‘i) = —3adF(p2 — Sa'^dpifi — a^dp^p (V.25) 

etc. 


This leads us to the following lemma: 
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Lemma V.2 




k=l 


(V.26) 


Proof. Having demonstrated it for the j = 1 case, we now proceed by indnction: 


= < 


L{(pj+i) = L{dt(pj) + L{d^0j) 

Ozt^j “I” -^dxz^j “1“ BOzz^j COz^j 
dt{L{(fj)) + dz + Adx<fj + Bdz<fj - CCpj) 

-{dzA)dx<^j - {dzB)dz<^j + {dzC)!^j 
9t(L((pj)) + dz{L{(:pj)) + aAzdx0j + aBzdz0j - aCtpj 
dt{L{^j)) + dz{L{^j)) 

+a(L((pj) - dt(pj - uodx0j - wodz^Pj) 
dt{L{(pj)) + dz{L{(pj)) + aL{(pj) - adp^Pj 

dt (- ELi {i)a’^dF<Pj-k) + dz (- Ei=i {i)a’^dF<Pj-k) 
{ +a (- Ei=i {l)a'^dF‘Pj-k) - adp^Pj 

.fco f a I O ,7? N (A ^k+l; 


= < 


= < 


(V.27a) 

(V.27b) 

(V.27c) 

(V.27d) 

(V.27e) 

(V.27f) 

(V.27g) 


“ + dz0j-k) - Y1 U (V.27h) 

fc—1 \ / /c^O \ / 


j 

-i: 

k=l 

^ O' + 1 


i+1 


dpipj+i-k Oi^dpipj+i-k 

a^dpipj+i.k - a^^'^dpip 


J J 

k [k-l 


E 

k=l 


k 


a 


'dF(pj+i_k - a^^^dpip 


L{0j+i) = 


E Pt ° 

k=i ' ^ 


(V.27i) 

(V.27j) 

(V.27k) 

(V.28) 


We can now nse this L{(p) to remove the normal derivative terms. Proceeding with 
each (pj in snccession, we hrst have 


Bdz<p = —dt0 — Adx(p + C(p + D . 


(V.29) 


Next, we have 


L{(pi) = dt0i + Adxipi + Bdzipi - C(pi = -adptp , 


(V.30) 
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which we combine with (V.29) to get 


Bdz0i = - Adxifi + C0i - adt(p - au^d^ip - awodz(p 

/ 

-dt(pi - Adx(pi + C(pi - a {I - woB~^) dt(p 
—a{uQl — WQB~^A)dx0—oiWQB~^{C(p + D^ . 


Lemma V.3 


Bd^ifj 


-dt(pj - Ad^ifj + C(pj - Ei=i {Cjoi'^Pk [(/ - wqB dt(pj-k 
< + {uqI - woB~^A) d^0j_k + woB~^C0j_k] 

—a^PjWoB~^D , 


for all j > 1, where the polynomial sequence Pj is defined recursively by 


Pj {woB-^) = ^ - E {woB-^) . 

Proof. As with the preceding lemma, having demonstrated the j = 1 
proceed by induction: 


= dt^fj+i + Adx^fj+i + BdzP>j+i — Cifj+i 




— Adx^j^i + Cipjj^i — ^ 


i+l + IN 


(y^dpfij+i-k 


t+i 

- E 

k=l 


k=l 


i 1 _ fc 


k 


a {dt + Uodx + Wodz) fij+i-k 


<h((Pj+i) - ELi {dtfij+i-k + uodxfij+i-k 

PwqB ^ {—dt^Pj+i-k — Adx<Pj+i-k + C(pj+i-k 

- ES“' [{I - WoB~^) dtfij+,_k-i 

+ {uqI - woB-^A) dxfij+i-k-i + woB-^Cfij+i_k-i] 

-a^+^-^Pj+i_kWoB~^D^'j - {dtfi + uodxfi 

+wqB~^ [-dtfi - Adxfi + Cfi + D^j 


(V.31) 

(V.32) 

(V.33) 
case, we 

(V.34) 

(V.35a) 

(V.35b) 

(V.35c) 
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' »(w+i) - Ei=i (':>* [(/ - »oS^') d,^i+,-k 

+ {uol - woB~^A) d^0j+i_k + woB-^C0j+i_k] 

Ei.i !(/ - 

+ {uqI - wqB^^A) dx0j+i-k-i + woB~^C0j+i-k-i] 

+woB-^ ELi - a^+i ((/ - woB-^) dt^p 

+ {uqI — wqB~^A) dxP + wqB~^ {C(p + 

^ $(^,+i) - ES [(/ - WoB-^) dt0,+,.k 

+ {uqI - wqB-^A) da,(fij+i-k + woB^^Cipj+i_k\ 

< +WoB-^ Ei=i ('t') [(/ - woB-^) dtPj+,iY^?Q) 

+ {uqI - woB~^A) dx0j+i-k-i + woB~^C0j+i-k-i] 


[ -aB^ (/ - ELi WoB-^D . 

Since the binomial coefficients are symmetric, we can replace Pj+i-k with Pk in the 
last line of the previons eqnation. The terms in the hnal gronp of parentheses are 
thns Pj+i by dehnition. As a temporary shorthand, let 


n 


(/ - wqB dtPj+i-k-i 

< 

+ {uqI - wqB-^A) d^pj+i-k-i + woB-^Opj+i-k-i 

aB^Pj+iw^B'^D . 


This gives ns 


Bdz^j-\-i \ 


<i>(^i+i) - Eis (':>'• [(/ - woB-') 

+ {uqI - woB~^A) d^(fj+i_k + woB~^C(pj+i_k] 

+i»oB-> Ei., E&y‘ (y‘) ('+i-‘)a‘*'-P'®(^j+i-‘-o 


(V.37) 


n 


We wonld like to replace the donble snmmation with one in terms of rather than 
aB’-, Let r = k + 1. Looking at the valnes assnmed hy k + I and /, and recognizing 
that I < k + I hy dehnition, we can rewrite this donble snmmation as 




r=ll=l 


I 
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A quick check shows that, not only are the indices mapped correctly, we also have 
the correct number of terms in the new summation: 


j j+l-k j+lr-l 

E E i = EEi = '-^ 

k=l 1=1 r=l1=1 ^ 


Also, we can simplify the binomial expansions by 


O' + 


Returning to our derivation, and reverting r back to /c, we have 


Bd^iPj+i = < 


= < 


= < 


- ES !(/ - 

+ {uol - wqB^^A) dx<fj+i-k + woB~^C(fj+i-k] 

[ +w„B-' eS Eti‘ (d‘) (?)a‘-Pie(ft+i-<=) - » 

[ eS Eti‘ (d‘) (?)«‘-Pie(ft+i-t) - Si 

Eili (/ - Eti‘ e(^j+i-i) 


V k 
i+1 




J + 1 \ ^ fc ; 


k=l 


k 


Bd^ipj+I = 


^t^j+i -^^x^j+i A Cipjj^i 

- Eit\ [(/ - woB-^) dt<fj+,_k 

+ {uol - wqB^^A) dx<fj+i-k + woB~^C(fj+i-k] 


-a- 


’+^Pj+iWoB-^D □ 


(V.38a) 

(V.38b) 

(V.38c) 

5d 

(V.38d) 

(V.39) 


We are now able to remove the normal derivative terms from the auxiliary variable 
equations. The auxiliary variable formulation for Fy is 


0j+i = {dt + 0j (V.40) 

(po = 0 

<Fj = 0 . 


105 



Using our formula for dz<fj, we have 


dt(fj + B ^ {-dt(pj - Ad^ipj + C 0 j 

- eLi (i)a*u l(^ - k'oB'’) 

+ (uqI - woB~^A) dx(pj-k + WQB~^C(pj-k] 
—a^Pj'WQB~^D^ . 


(V.41) 


Collecting all dt terms on the left-hand side with the (pj+i term, leaving all the other 

terms on the right-hand side, and left-multiplying everything by 5, we have 

/ 

/ ■\ — Adx^-Pj C(pj 

Y:l-i(i)a'"PkU-woB-Adt0.k -A. 

^ = -j:i=AiyPk[{uoi-woB-^A)dx0,.k 

+ {I-B) dt 0 i + B 0.+1 ^ 

+woB-^C0j-k] - a^PjWoB~W . 

(V.42) 

In characteristic form, we left-multiply both sides by R~^ and insert RR~^ as needed 


to get 

Ei.i 

+ (/ — A) dtCj + A^+i 

for j G 1, 2,..., J — 1, where 


—Adxij + Cij 

- ELi i^^o.^Pk [{uol - woA-^A) dxij-k 
+woA~^C^j-k - a^PjWoA~^D , 

(V.43) 
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A = R-^AR 

P, = P-^P,P (V.44) 

and ^j, C, and D are defined as in (V.9). As with the simple case, we nse natnral 
bonndary conditions for the ontgoing characteristics and the trnncation condition for 
the incoming. Hence, onr anxiliary variable system is 



(/-A)a,6 + Aei 


-Ad^io - Ad^io + {A{d,R-^)R 
~Adx^o + + D 


E 


j 

k=l 


f -AdxCj + CC, 


g)aVa/-woA-i)aV_fc 1 

> = < 

-EJ.1 

+ (/ — A) dtij + A^-+i j 


-WqA-^ A)dxij-k 



^ +woA ^Cij-k 

- 

[e 

/ = 0 

1 

5 



{uqI 


a^PjWoA 

(V.45) 


where we again nse the + and — snbscripts to denote ontgoing and incoming char¬ 
acteristics, respectively. Note also that, in order to have the term in the first 
eqnation, we nse the prodnct rnle and the ^-dependence of R~^: 

dz(o = dz (P^ Vo) = R~^d^(pQ -F {d^R~^^ (po . 

In the absence of gravity, a = 0, and (V.45) is analogons to (V.IO). 

We lose some compntational efficiency with this implementation. Eq. (V.IO) 
defined a bi-diagonal matrix system, which reqnires 0{J) operations per time step 
with an efficient matrix solver. On the other hand, (V.45) is lower-Hessenberg and 
reqnires 0{J^) operations per time step. 

When we extend this derivation to the open domains defined in the previons 
section, we get the following system nsing N and T instead of A and B (compare to 
(V.19)): 

-AaV- fdxi + (A {dnR-^) R 

+f{drR-^)R + C)i + D\^ 
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{I-A)dtio + Hi = -f{dM-{drR-^)Rio)+C(o + D 

[ -f{d^(j-{d^R-^)RQ+C(j 


ELi (/ - noA-i) dtCj-k 
+ (/ — A) dtCj + A0j+i 


> = < 


0 


= 0 


-Ei.ig)a''-P4(T-„/ 

) (Brij-t - (drR-') Rij-k) 

[ +noA~^C^j-k - &PjnoA~^D 

(V.46) 


where 


Pk = Pk{noN-^) 


A — diag no + Co no no no - co 
N = R-^AR 


f = R-^TR 


a 


a = 


V2 ■ 


In addition, we must replace po in the eigenvector matrices R and R~^ with the 
^-dependent p. 

About the Polynomial Sequence. A quick literature search [2, 60, 103, 118] 
revealed no matches to the polynomial sequence 


n—1 


Pn{x) = 1 - Xj • 


The hrst hve polynomials of this sequence are 


Pi = 

1 


P2 = 

1 -2a: 


Ps = 

1 — 6a: + 6a:^ 


Pa = 

1 — 14a: + 36a:^ - 

-24a:^ 

n = 

1 - 30a: + ISOa:^ 

- 240a:2 + 120a:^ 


The hrst three are constant multiples of the Bernoulli polynomials [103], but not the 
subsequent ones. We have not found any other matches to this sequence. However, 
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Dr. Pante Stanica at the Naval Postgraduate School observed [104] that the coeffi¬ 
cients of Pn can be generated directly, without use of the recursion sequence, by the 
formula 

where |^| denotes the Stirling number of the second kind [46]. Other properties of 
these polynomials, such as orthogonality and a generating function, may be explored 
in future research. 

2. Numerical Examples 

Let us consider a few numerical examples, using our standard simulation. 
First, we look briefly at the system in a semi-inhnite channel, open on top (see 
Fig. 6), subject to gravitational forces. There is no real advection, but we dehne our 
false boundary advection wq = = 1.0417^ to keep the matrices non-singular. 

Table XXII shows the error norms resulting from the cases J G 1... 10. 



J 

E, 

Eu 

Ev 

Ep 

1 

0.073179 

0.08778 

0.11305 

0.049617 

2 

0.073086 

0.087435 

0.11249 

0.04944 

3 

0.07298 

0.08704 

0.11197 

0.049261 

4 

0.072874 

0.086646 

0.11146 

0.049083 

5 

0.07277 

0.086255 

0.11095 

0.048907 

6 

0.072667 

0.085866 

0.11045 

0.048732 

7 

0.072565 

0.08548 

0.10995 

0.048558 

8 

0.072464 

0.085096 

0.10945 

0.048386 

9 

0.072364 

0.084716 

0.10896 

0.048215 

10 

0.072264 

0.084337 

0.10848 

0.048046 

Som. 

0.092721 

0.10278 

0.22495 

0.092238 


Table XXII. Error norms (IV. 14) for Givoli-Neta NRBCs, J G 1... 10, gravity- 
influenced system (IV.31), semi-inhnite vertical channel, no advection. The error 
norm from using a Sommerfeld boundary condition is included for comparison. 


Next, we consider the open-air domain (Fig. 22). We begin hrst with “zero” 
advection (false boundary advection uq = = 1.3889^, wq = = 1.0417^, 
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chosen so that uq 7 ^ Wq to keep N non-singular). Fig. 39 shows the state variable p 
at the end of the run for J = 10. Table XXIII shows the error norms (IV. 14) for each 


Pressure perturtation, reference solution 



Pressure perturtialion 
100 
60 
60 
40 
20 



I 

i 


Tnjncated reference solution 

100 I- 


|50 80 

60 
40 

I 20 

•50 





Figure 39. The solution for the pressure perturbation p using J = 10, Givoli-Neta 
NRBCs in an open half-plane subject to gravitational forces with no advection. 
(Top) Reference solution; the area corresponding to the computed solution is con¬ 
tained in the bottom-center box. (BL) Computed solution using NRBCs. (BC) Ref¬ 
erence solution domain corresponding to NRBC solution domain. (BR) Delta between 
computed and reference solutions, with error norm computed by (IV. 14). 


state variable as J goes from 1 to 10. 

For our third example, we use the non-zero mean flow Uq = 100^ (maintaining 
the same “false” Wq on the boundary), we get the error norms in Table XXIV. 

In both cases, we observe that there is very little improvement when we in¬ 
crease J. This property of auxiliary variable NRBCs has been noted before in other 
contexts [36] and is the result of discretization errors induced by converting the nor¬ 
mal derivative terms to tangential. Higher-order discretization schemes may improve 
these results [55], but for the purpose of this dissertation, it is sufficient to demon¬ 
strate how to use the auxiliary variable NRBC methods with the linearized Euler 
equations. Furthermore, we note by the bottom row of Tables XXII-XXIV that the 
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J 

E, 

Eu 

Ey 

Ep 

1 

0.30742 

0.48811 

0.20486 

0.21176 

2 

0.30615 

0.48976 

0.20039 

0.20991 

3 

0.30532 

0.49034 

0.19818 

0.2087 

4 

0.30472 

0.49024 

0.19688 

0.20784 

5 

0.30426 

0.48972 

0.19598 

0.20717 

6 

0.30387 

0.48896 

0.19524 

0.20662 

7 

0.30353 

0.48805 

0.19459 

0.20613 

8 

0.30321 

0.48706 

0.19399 

0.20568 

9 

0.30291 

0.48602 

0.19342 

0.20526 

10 

0.30261 

0.48497 

0.19286 

0.20485 

Som. 

0.30224 

0.85386 

0.30744 

0.29872 


Table XXIII. Error norms (IV. 14) for Givoli-Neta NRBCs, J G 1... 10, gravity- 
influenced system (IV.31), open-air domain, no advection. The error norm from 
using a Sommerfeld boundary condition is included for comparison. 

auxiliary variable method is signihcantly more accurate than the basic Sommerfeld 
(hrst-order Higdon) condition. 

Note: Examination of our characteristic systems (V.45) and (V.19) reveals 
that, in a uniform domain (zero perturbations in all state variables), the presence 
of the non-zero D on the right-hand side will generate non-zero results at the next 
time step. Physically, this should not happen; it is a consequence of replacing the 
normal derivatives with tangential derivatives. To avoid this problem, a check should 
be made at each point on the boundary prior to applying the boundary condition: If 
the solution at the point, its immediate neighbors, and the corresponding auxiliary 
variables are all zero, then the solution at the point should remain zero. 

D. SUMMARY 

In this chapter we have modihed the Givoli-Neta auxiliary variable method for 
implementation in a hrst-order system. Although our original context is the linearized 
Euler equations, the linear matrix equation can be extended easily to any hrst-order 
system. The gravity-inhuenced derivation is specihc to the linearized Euler equations. 
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J 

E, 

Eu 

Ey 

Ep 

1 

0.20976 

0.25767 

0.24514 

0.18111 

2 

0.20923 

0.25683 

0.24286 

0.18049 

3 

0.20866 

0.256 

0.24061 

0.17987 

4 

0.20813 

0.25517 

0.2384 

0.17925 

5 

0.2076 

0.25435 

0.23623 

0.17863 

6 

0.20703 

0.25353 

0.2341 

0.17802 

7 

0.20647 

0.25272 

0.232 

0.17741 

8 

0.20596 

0.25191 

0.22994 

0.17681 

9 

0.20546 

0.25112 

0.22793 

0.17621 

10 

0.20492 

0.25036 

0.22596 

0.17562 

Som. 

0.25189 

0.34105 

0.27768 

0.24917 


Table XXIV. Error norms (IV. 14) for Givoli-Neta NRBCs, J G 1... 10, gravity- 
influenced system (IV.31), open-air domain, left-to-right advection. The error norm 
from using a Sommerfeld boundary condition is included for comparison. 


but the rest of the implementation is fully portable to any other linear hrst-order 
system. 

In the next chapter we consider a third class of NRBC, one which was recently 
developed in the literature, which we shall extend to our equation system. 


112 



VI. HAGSTROM-WARBURTON NRBCS 
FOR THE LINEARIZED EULER EQUATIONS 


A. INITIAL IMPLEMENTATION FOR FIRST-ORDER 
SYSTEMS 
1. Derivation 

We now consider a third class of NRBC, the Hagstrom-Warbnrton anxiliary 
variable techniqne, and apply it to the linearized Enler eqnations. We proceed simi¬ 
larly to Sec. 6 of [57], which docnmented the concept’s extension to symmetric hrst- 
order systems. While onr system can be made symmetric [117], it is actnally not 
necessary to do so, and some nnmerical experiments have shown that the resnlts are 
either the same as the following straightforward approach, worse, or nnstable, depend¬ 
ing on the domain conhgnration. The vector form of (III.33) is a fairly straightforward 
extension. 


dt(pi = aodt<f + dfi(p 

ajdt(pj+i - dfi0j+i = ajdt(pj+i + 8^0j (VI. 1) 

0J+1 = 0 

We remove the wave speed c from the eqnations becanse the wave speeds will instead 
be replaced by the eigenvalnes of the coefficient matrix corresponding to the normal 
derivative. As with the preceding chapter for the Givoli-Neta formnlation, we begin 
with a nsefnl lemma: 

Lemma VI.1 Let L{(p) = dt(p + Adx(p + Bdytf — C(p. The auxiliary variables ipj 
satisfy 

m) = 0 (VI.2) 

for all j > 1. The only exception is when both of the following conditions apply: 

• n = ±/c 

• ^ 0 
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Proof. Apply the differential operator L to (VI.la), giving 


L [dtipi) = L {aodt0 + dfi0) 

= (aodt(L{D))+dn{L{D))) 

dt{L{ip,)) = 0 (VI.3) 

Since 0i is initially zero, L{0i) is also zero, and thns it is always zero. Proceeding by 
indnction we now apply L to (VI. lb) and get 

L {oijd^ipdfi^p L {o,jdi(pj 

= ajdt{L{(pj)) + dri{L{(pj)) 

= 0 (VI.4) 

We then integrate along the incoming characteristic to get L{(pjj^i), which mnst be 
zero becanse the incoming characteristics are zero □ 

Once again, as with Lemma V.l, the proof relies on the fact that L{dfi(p) = 
dfi [L {(p)). On the top and bottom of a gravity-inffnenced domain, that is not trne, 
becanse the coefficient matrices also depend on We will address this issne in Sec. C. 

As with the Givoli-Neta formnlation, we wonld like to get rid of the normal 
derivative terms in (VI. 1). For definiteness, the following derivation is applied to a 
right-hand-side bonndary, that is, n = i. Simply make the appropriate changes for 
the other bonndaries (except a vertical bonndary in the presence of gravity, as we 
have noted). Left-mnltiply (VI.la) by A and snbtract 0 = L{(p) — D to get 

Adt(pi = (ooA — I) dt(p — Bdyip + Cip -f D (VI.5) 

Similarly, left-mnltipIy (VI.lb) by A and add L{(pj^i) = —L{pj) = 0 to get 

{ajA + I) dtPj+i + BdyPj+i - Cpj+i = {ajA - I) dtPj - BdyPj + Cpj (VI.6) 

Consider the eigenvalne and eigenvector matrices for the Givoli-Neta NRBCs 
(V.7). Left-mnitiply (VI.5) and (VI.6) by R~^ and make the snbstitntion = R~^(pj. 
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We get 


— (aoA — I) dt^o — Bdy^o + C^o + B 

(a,A + I) dt{j+i + - Cg+i = (a,A - I) dtC, - Bdyi, + Cg , (VI.7) 

where B, C and D are all defined by (V.9). We also insert R~^R as needed to 
change the (pj terms to ^j’s. We have a self-contained system of eqnations dehned 
solely on the bonndary. Hence, with the bonndary pertnrbations initially zero, these 
eqnations will always remain zero. So we need to incorporate interior valnes into 
onr formnlation in order to have something more than a rather convolnted Dirichlet 
condition. Following the approach in [57], we replace the fonr-eqnation trnncation 
condition (pj+i = 0 with a more characteristic-based approach. For the ontgoing 
characteristics, we nse the interior scheme 

'dtio + Aa^go + Bdyio = Cio + d]^ , (VI.8) 

where the -f snbscript indicates that we only consider the resnlting eqnations corre¬ 
sponding to the ontgoing characteristics. This gives ns one or three eqnations (de¬ 
pending on the sign of Uq) for onr system of nnknowns. For the remaining nnknowns, 
we apply the trnncation condition to the incoming characteristics only, i.e., 

[ 6 h-i = 0]_ (VI.9) 

These two eqnations, in conjnnction with (VI.7), dehne a system of 4( J + 2) eqnations 
for 4( J -|- 2) nnknowns ^j. Collecting the time derivative terms onto one side gives ns 
the explicit system: 


dt^o 

A^ig + (/ — agA) dtio 
(/ + cijA) dt^j+i + (/ — cijA) dt^j 

[gj+i 


-AdJo - Bdyjo + {b {d,R-^) R + C)^o + d]^ 
-Bdyjo + {b (d,R-^) R + C)io + D 
-Bdy\J+i + (b {d,Rj R + C) g+i 
^ -Bdy\J + {B{d,RjR + C)iy 
0]_ , (VI.IO) 
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where the subscript y\z denotes a partial derivative with respect to either y or de¬ 
pending on the domain. The dzR~^ term is non-zero only when gravity is considered. 

Note that the presence of A on the left-hand side of (VI. 10b) requires that we 
again have A (and thus A) non-singular, as with the Givoli-Neta NRBC. We solve the 
problem here the same way. Impose a false advection uo = e ^ 0 small enough that 
the system will not drift across a single grid point over the duration of the simulation; 
thus, 

St 

ll«oll < 5 — (VI.ll) 

^^max 

where Sx is the grid spacing, tmax is the duration of the simulation, and we include a 
factor of two to ensure that there will be no drift even with round-off errors. 

2. Numerical Examples 

Let us consider a few numerical examples. We use our usual inhnite-channel 
example, considering the same six cases as Sec. V.3, always with J G 1... 10. Ta¬ 
ble XXV shows the state variable error norms for the basic system with no advection. 
Table XXVI shows the error norms for the basic system with left-to-right advec¬ 
tion. Table XXVII shows the error norms for the Coriolis-influenced system with 
no advection. Table XXVIII shows the error norms for the Coriolis-influenced sys¬ 
tem with left-to-right {uq = 100^) advection. Table XXIX shows the error norms 
for the gravity-influenced system with no advection (recall that we can handle grav¬ 
ity when our open boundaries are on the left and/or right, but not on the top or 
bottom). Table XXX shows the error norms for the gravity-influenced system with 
left-to-right advection. For comparison, each table includes the error norm resulting 
from using a Sommerfeld condition. For the “no advection” cases we use the false 
advection uq = 1.3889^. In every case, we make the simplihcation aj = 1 for all j. 
Hagstrom et ah [53], citing an analysis by Diaz and Joly [23], show that this choice is 
always adequate. We will also make this simplihcation in all subsequent derivations 
for these NRBCs. Fig. 40 shows the state variable p at the end of the run for J = 10 
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for the gravity-influenced system with left-to-right advection. Note again the strong 
reflections which accompany the characteristic-based implementations [72], 


Pressure pefturt>atior^. reference solution 



-50 0 50 100 150 200 


Pressure perturbation Truncated reference solution Solution delta 



SO 100 SO 100 SO 100 


Error nomv-0.12d2 


Figure 40. The solution for the pressure p using J = 10, Hagstrom-Warburton 
NRBCs, gravity-influenced system (IV.31), inhnite horizontal channel, left-to-right 
advection. (Top) Reference solution; the area corresponding to the computed so¬ 
lution is contained between the vertical lines. (BL) Computed solution using NR¬ 
BCs. (BC) Reference solution domain corresponding to NRBC solution domain. 
(BR) Delta between computed and reference solutions, with error norm computed by 
(IV.14). 
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J 

E, 

Eu 

Ey 

Ep 

1 

0.12441 

0.25176 

0.10978 

0.12441 

2 

0.12353 

0.24955 

0.10744 

0.12353 

3 

0.12264 

0.24735 

0.10619 

0.12264 

4 

0.12175 

0.24519 

0.10511 

0.12175 

5 

0.12086 

0.24305 

0.10407 

0.12086 

6 

0.11999 

0.24095 

0.10304 

0.11999 

7 

0.11912 

0.23887 

0.10203 

0.11912 

8 

0.11827 

0.23683 

0.10103 

0.11827 

9 

0.11742 

0.23482 

0.10004 

0.11742 

10 

0.11658 

0.23283 

0.099065 

0.11658 

Som. 

0.21847 

0.53616 

0.12503 

0.21847 


Table XXV. Error norms (IV. 14) for Hagstrom-Warburton NRBCs for basic system 
(IV. 1) with J G 1... 10 in an in fini te channel with no advection. Error norms from 
nsing a Sommerfeld radiation condition are inclnded for comparison 


J 

E, 

Eu 

Eu 

Ep 

1 

0.13788 

0.2047 

0.14207 

0.13788 

2 

0.1367 

0.20297 

0.13915 

0.1367 

3 

0.13553 

0.20126 

0.13632 

0.13553 

4 

0.13437 

0.19958 

0.13358 

0.13437 

5 

0.13322 

0.19791 

0.13092 

0.13322 

6 

0.13209 

0.19626 

0.12835 

0.13209 

7 

0.13096 

0.19462 

0.12586 

0.13096 

8 

0.12985 

0.19301 

0.12344 

0.12985 

9 

0.12874 

0.19141 

0.1211 

0.12874 

10 

0.12765 

0.18983 

0.11884 

0.12765 

Som. 

0.21817 

0.29414 

0.17727 

0.21817 


Table XXVI. Error norms (IV. 14) for Hagstrom-Warbnrton NRBCs for basic system 
(IV. 1) with J G 1... 10 in an in fini te channel with left-to-right advection. Error 
norms from nsing a Sommerfeld radiation condition are inclnded for comparison 
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J 

E, 

Eu 

Ev 

Ep 

1 

0.12441 

0.25176 

0.10979 

0.12441 

2 

0.12353 

0.24955 

0.10745 

0.12353 

3 

0.12264 

0.24736 

0.10621 

0.12264 

4 

0.12175 

0.24519 

0.10513 

0.12175 

5 

0.12086 

0.24306 

0.10408 

0.12086 

6 

0.11999 

0.24095 

0.10306 

0.11999 

7 

0.11912 

0.23888 

0.10204 

0.11912 

8 

0.11827 

0.23684 

0.10104 

0.11827 

9 

0.11742 

0.23482 

0.10005 

0.11742 

10 

0.11658 

0.23284 

0.09908 

0.11658 

Som. 

0.21847 

0.53616 

0.12503 

0.21847 


Table XXVII. Error norms (IV. 14) for Hagstrom-Warburton NRBCs for Coriolis- 
influenced system (IV. 15) with J G 1... 10 in an infinite channel with no advection. 
Error norms from using a Sommerfeld radiation condition are included for comparison 


J 

E, 

Eu 

Eu 

Ep 

1 

0.13855 

0.31144 

0.21862 

0.13855 

2 

0.13753 

0.30948 

0.21789 

0.13753 

3 

0.13654 

0.3076 

0.21725 

0.13653 

4 

0.13557 

0.3058 

0.21668 

0.13557 

5 

0.13464 

0.30407 

0.21617 

0.13464 

6 

0.13372 

0.3024 

0.21573 

0.13372 

7 

0.13284 

0.3008 

0.21534 

0.13284 

8 

0.13198 

0.29926 

0.215 

0.13197 

9 

0.13113 

0.29777 

0.21471 

0.13113 

10 

0.13032 

0.29633 

0.21446 

0.13032 

Som. 

0.1929 

0.32104 

0.15953 

0.1929 


Table XXVIII. Error norms (IV. 14) for Hagstrom-Warburton NRBCs for Coriolis- 
influenced system (IV. 15) with J G 1... 10 in an inhnite channel with left-to-right 
advection. Error norms from using a Sommerfeld radiation condition are included for 
comparison 
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J 

E, 

Eu 

Ev 

Ep 

1 

0.12896 

0.25097 

0.11173 

0.12835 

2 

0.12807 

0.24879 

0.10934 

0.12745 

3 

0.12714 

0.24663 

0.10807 

0.12652 

4 

0.12622 

0.24449 

0.10697 

0.1256 

5 

0.1253 

0.24238 

0.10591 

0.12469 

6 

0.12439 

0.24031 

0.10487 

0.12379 

7 

0.12349 

0.23826 

0.10384 

0.12289 

8 

0.12259 

0.23624 

0.10283 

0.12201 

9 

0.12171 

0.23425 

0.10183 

0.12113 

10 

0.12084 

0.2323 

0.10084 

0.12026 

Som. 

0.22442 

0.52603 

0.12775 

0.22374 


Table XXIX. Error norms (IV. 14) for Hagstrom-Warburton NRBCs for gravity- 
influenced system (IV.31) with J G 1... 10 in an infinite channel with no advection. 
Error norms from using a Sommerfeld radiation condition are included for comparison 


J 

E, 

Eu 

Ey 

Ep 

1 

0.13899 

0.20395 

0.14268 

0.13852 

2 

0.13779 

0.20223 

0.13977 

0.13733 

3 

0.1366 

0.20053 

0.13696 

0.13615 

4 

0.13543 

0.19885 

0.13423 

0.13498 

5 

0.13426 

0.19718 

0.13158 

0.13382 

6 

0.13311 

0.19554 

0.12902 

0.13267 

7 

0.13197 

0.19391 

0.12654 

0.13154 

8 

0.13084 

0.1923 

0.12414 

0.13041 

9 

0.12972 

0.19071 

0.12181 

0.1293 

10 

0.12861 

0.18913 

0.11956 

0.1282 

Som. 

0.21869 

0.28969 

0.17855 

0.21802 


Table XXX. Error norms (IV. 14) for for Hagstrom-Warburton NRBCs for gravity- 
influenced system (IV.31) with J G 1... 10 in an inhnite channel with left-to-right 
advection. Error norms from using a Sommerfeld radiation condition are included for 
comparison 
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B. CORNER CONDITIONS 
1. Derivation 

We can use the same methods as in Sec. V.l to apply these NRBCs to the 
corners of an open domain, in the absence of gravity. Using the same notation as 
before, our characteristic-based NRBC scheme on three open sides of a half-plane is 
given by 


dtC+ = -AdnC+ - Tdr 0 + 
Adti, + {I-A)dti = -fdr0 
{I + A)dt(j^, + {I-A)dti, = 

ij+i- = 0 


where we use 


dn 

Left : < dr 

T 

.. 

dn 

Right : < dr 
T 

.. 

dn 

Top : < dr 

T 

dn 

Top-left : < dr 

T 

dn 

Top-right : dr 

T 


-dx 

dz 

B 

dx 

dz 

B 

dz 

dx 

A 


(—dx + dz) /-\/2 
(dx + dz) / 

(A + B)/V2 

(dx + dz) / \f2 

(dx - dz) /v^ 

(/I - B) /y/2 


(VL12) 


(VL13) 
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2. Numerical Examples 

Let us now consider some examples. Since we cannot include the effects of 
gravity on the top of the domain, we consider only four cases: The basic system 
and the system influenced by Coriolis forces, with and without horizontal advec- 
tion. We use the same false advection as in Sec. 2 (horizontal) and V.2 (vertical). 
Table XXXI shows the error norms for J G 1... 10 for the basic system with no ad¬ 
vection. Table XXXII shows the error norms for J G 1... 10 for the basic system with 
left-to-right advection. Table XXXIII shows the error norms for J G 1... 10 for the 
Coriolis-influenced system with no advection. Table XXXIV shows the error norms 
for J G 1... 10 for the Coriolis-influenced system with left-to-right advection. Fig. 41 
shows an example of the horizontal velocity perturbation u for the basic no-advection 
test. Note the reflections are stronger at the bottom corners than the top. These 
reflections are the result of waves reflecting off the hard wall on the bottom before 
impacting the left and right boundaries; thus, their incidence angles are greater than 
45°. At the top corners, these waves pass through the open top with the smaller inci¬ 
dence angle, so there is less reflection evident. If we plot the values of Table XXXIV, 
we see that the error norms decrease exponentially as J increases by 2 (Fig. 42, but 
note that v increases with J > 7). The reason for the saw-tooth pattern is unknown, 
but it is only present in the Coriolis-l-advection case; the other cases show exponential 
reductions with each increase in J. 
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-50 0 SO 100 ISO 200 

X-direction velocity perturbation Truncated reference solution Solution delta 



Figure 41. The solution for the horizontal velocity u using J = 10, Hagstrom- 
Warburton NRBCs, basic system (IV.1), open half-plane, no advection. (Top) Ref¬ 
erence solution; the area corresponding to the computed solution is contained in the 
bottom-center box. (BL) Computed solution using NRBCs. (BC) Reference solution 
domain corresponding to NRBC solution domain. (BR) Delta between computed and 
reference solutions, with error norm computed by (IV. 14). 


J 

E, 

Eu 

Ey 

Ep 

1 

0.17603 

0.30011 

0.13089 

0.17604 

2 

0.17485 

0.29724 

0.12848 

0.17485 

3 

0.17364 

0.29438 

0.12699 

0.17365 

4 

0.17244 

0.29157 

0.12565 

0.17245 

5 

0.17125 

0.28881 

0.12435 

0.17126 

6 

0.17007 

0.2861 

0.12307 

0.17007 

7 

0.1689 

0.28345 

0.12182 

0.1689 

8 

0.16773 

0.28085 

0.12059 

0.16774 

9 

0.16658 

0.2783 

0.11938 

0.16658 

10 

0.16543 

0.2758 

0.11819 

0.16543 

Som. 

0.33324 

0.71942 

0.2414 

0.33325 


Table XXXI. Error norms (IV. 14) for Hagstrom-Warburton NRBCs for basic system 
(IV. 1) with J G 1... 10 in an open half-plane with no advection. Error norms from 
using a Sommerfeld radiation condition are included for comparison 
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J 

E, 

Eu 

Ey 

Ep 

1 

0.16823 

0.20592 

0.18956 

0.16823 

2 

0.16661 

0.20403 

0.18543 

0.16661 

3 

0.16501 

0.20217 

0.18143 

0.16501 

4 

0.16342 

0.20034 

0.17755 

0.16342 

5 

0.16186 

0.19854 

0.1738 

0.16186 

6 

0.16031 

0.19675 

0.17016 

0.16031 

7 

0.15878 

0.19499 

0.16664 

0.15878 

8 

0.15727 

0.19326 

0.16323 

0.15727 

9 

0.15577 

0.19154 

0.15993 

0.15577 

10 

0.15429 

0.18985 

0.15673 

0.1543 

Som. 

0.2675 

0.32527 

0.2411 

0.2675 


Table XXXII. Error norms (IV. 14) for Hagstrom-Warburton NRBCs for basic system 
(IV. 1) with J G 1... 10 in an open half-plane with left-to-right advection. Error 
norms from nsing a Sommerfeld radiation condition are inclnded for comparison 


J 

E, 

Eu 

Eu 

Ep 

1 

0.17604 

0.30024 

0.13095 

0.17604 

2 

0.17484 

0.29737 

0.12852 

0.17484 

3 

0.17364 

0.29451 

0.12704 

0.17365 

4 

0.17243 

0.29169 

0.12569 

0.17243 

5 

0.17125 

0.28893 

0.12439 

0.17126 

6 

0.17006 

0.28622 

0.12311 

0.17006 

7 

0.16889 

0.28356 

0.12186 

0.1689 

8 

0.16772 

0.28096 

0.12063 

0.16772 

9 

0.16657 

0.2784 

0.11943 

0.16658 

10 

0.16541 

0.27591 

0.11824 

0.16542 

Som. 

0.33324 

0.71943 

0.2414 

0.33325 


Table XXXIII. Error norms (IV. 14) for Hagstrom-Warbnrton NRBCs for Coriolis- 
inflnenced system (IV. 15) with J G 1... 10 in an open half-plane with no advection. 
Error norms from nsing a Sommerfeld radiation condition are inclnded for comparison 
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J 

E, 

Eu 

Ey 

Ep 

1 

0.20445 

0.35374 

0.261 

0.20444 

2 

0.20472 

0.35221 

0.26198 

0.20471 

3 

0.20277 

0.34939 

0.26073 

0.20276 

4 

0.20307 

0.34806 

0.26177 

0.20306 

5 

0.20117 

0.34542 

0.26059 

0.20116 

6 

0.20152 

0.34424 

0.26168 

0.2015 

7 

0.19966 

0.34175 

0.26057 

0.19965 

8 

0.20003 

0.34072 

0.26169 

0.20002 

9 

0.19822 

0.33836 

0.26064 

0.19821 

10 

0.19862 

0.33744 

0.26179 

0.19861 

Som. 

0.26015 

0.34856 

0.20456 

0.26014 


Table XXXIV. Error norms (IV. 14) for Hagstrom-Warburton NRBCs for Coriolis- 
influenced system (IV. 15) with J G 1... 10 in an open half-plane with left-to-right 
advection. Error norms from using a Sommerfeld radiation condition are included for 
comparison 


Density perturbation error norm X-velocity perturbation error norm 



Y-velocity perturbation error norm Pressure perturbation error norm 




Figure 42. Logarithmic plot of error norms (IV. 14) for Hagstrom-Warburton NRBC, 
J G 2 ... 10, open half-plane, Coriolis, left-to-right advection. (TL) Error norm for p. 
(TR) Error norm for u. (BL) Error norm for v. (BR) Error norm for p. 
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C. INCOMPATIBILITY WITH GRAVITY 


Problems arise, however, when we try to incorporate the effects of gravity on 
the top bonndary. (We can inclnde it on P/, and Pr withont any difficnlty.) Using 
the same approach as for the Givoli-Neta NRBCs, we apply the differential operator 
L to both sides of (VI.la), beginning the process of removing the normal derivatives: 


L (5t(pi) = L (ao5t(p + 






aodt{L{(p)) + L{d^(p) 

aodtD + dt^(f + - Cd^(p 

/ 

dz {dt(p + Ad^ip + Bdz0 - C(p) 

- {dzA) (p - {dzB) (p + {dzC) (p 
d, (L m - {3^)0- {d,B)0+ {d,C) (p . (VI. 14) 


Even with the simplifying assnmption of an exponentially-decaying atmosphere, we 
still have 


dt{,L [tpi)) 


dt{L {0i)) 


dzD + aAzip + aBz0 — aC(p 

-aD + a (L {tp) - dt0 - Uod^0 - Wodz<p) 

-aD + aD - a {dt(p -f uodx(p + wodz(p) 

-adF(p . (VI. 15) 


To continne this derivation, we mnst integrate the Lagrangian flow derivative of (p over 
time. The resnlt is then applied to find L (<^ 2 ), with increasingly convolnted combina¬ 
tions of Lagrangian derivatives and time integrals. This approach is not satisfactory. 
Hence, we reqnire a different approach if we wish to incorporate gravitational effects 
into the Hagstrom-Warbnrton NRBCs. We have not fonnd snch an approach, bnt 
fntnre research may reveal a snitable techniqne. 
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D. WAVE-LIKE IMPLEMENTATION 
1. Derivation 

The results from this implementation are not bad, but the improvement with 
larger J is slight. Furthermore, the reflections downstream of the advection are trou¬ 
blesome. Perhaps we can do better. Looking at (VI. 1), we see a hint of the problem. 
Although Hagstrom and Warburton do not provide a physical interpretation of their 
auxiliary variable formulation, inspection of (VI.1) shows that we can conceive of the 
second equation in the following manner: The outgoing characteristic of (pj is in some 
sense paired with the incoming characteristic of Pj+i- However, when we use these 
characteristic variables for a hrst-order system, we contradict this interpretation of 
the auxiliary variables, since each variable has only one characteristic (either incom¬ 
ing or outgoing), not two. If we could contrive a second characteristic for each state 
variable, it might improve our results. Let us instead apply the NRBC to each state 
variable individually. Instead of Lemma VI.l, we use the following: 


Lemma VI.2 Each state variable p G {p, has a solution which satisfies the 

acoustic wave equation 


daP = , 


(VL16) 


where 


Co + Uo 

on Te 

Co — Uq 

onTw 

Co + Vo 

on Tat 

o 

1 

o 

on T^ 


(VI. 17) 


The derivation is given in Appendix C. With this fact, we can then implement the 


auxiliary variables exactly as given in [57], taking each state variable individually. 
Upon removing the normal derivatives, we are left with the following system of equa¬ 


tions: 


—a^dtp + dtPi 
2ai(l - al)dttp 
-|-(1 -[- -[- 2aQai)dttPi ' 

V(1 ~ (^\)dttP2 


= CpdxP 

/ 

cl{2aidyyP 
+ dyyPl + dyyP2) 
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aj{l- 

-\-{aj-i + cij)[l + aj-iaj)dtt‘fj 
+aj_i(l — aj)du‘fj+i 

ap{l- al,_^)dtt(pp~i 
+ {ap-i + cip){l + ap-iap)dtt<fp 


> = 


> = 


cl{aAy^j-i 

+ (aj-i + aj)dyy(pj 
+aj-idyy(pjpi) 

cliapdyyifp-i 
+ (ap_i + ap)dyy(pp) 


(VL18) 


where = Co + Uq. When we make the simplihcation a, = 1, this system reduces to 


-dtif + dtifi 
^du(pi 

4^dtt(Pj 

^dtt<fp 


Cpdx<f 

/ 

Cl{‘^dyyip 

< 

+ dyy(fl + dyy(f2) 

/ 

+2dyy(pj 
+ dyy(pjpl) 

\ 

/ 

cl{dyy^P-l 

< 

-\-2dyy(pp) , 


(VL19) 


This system is the NRBC on Tp. On Tw replace dx^p in (VI.19a) with —dx^p, and 
use Cxu = Co — Uq- 


2. Numerical Examples 

If we run the same infini te-cha nn el example as in Sec. 1.2, using the wave-based 
auxiliary variable approach from the preceding section, we get the results illustrated in 
Fig. 43 with the error norms for all state variables given in Table XXXV. If we use the 
non-zero mean flow uq = 100^, we get the error norms in Table XXXVI. Comparing 
Tables XXXV and XXXVI to Tables XXV and XXVI, respectively, we see that this 
new version’s error norms are approximately half the old version’s. We note also that 
the new method’s error norms show almost no improvement for J > 5. This “error 
floor” has been observed in other auxiliary variable implementations for scalar wave 
equations [36]. Although the old method shows continual (albeit slow) improvement, 
the error norms do not match those of the new method until J ~ 500, which runs far 
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Pressure pefturt>a(ion Truncated reference solution Solution dela 
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60 
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Error norm-0-059238 



0 


I 


Figure 43. The solution for the pressure p using J = 10, wave-like Hagstrom- 
Warburton NRBCs in an inhnite channel, basic system with no advection. (Top) Ref¬ 
erence solution; the area corresponding to the computed solution is contained between 
the vertical lines. (BL) Computed solution using NRBCs. (BC) Reference solution 
domain corresponding to NRBC solution domain. (BR) Delta between computed and 
reference solutions, with error norm computed by (IV. 14). 

slower than the new method’s J = 5. While this result is promising, experiments with 
longer integration times show that this wave-like implementation is less stable than 
the characteristic-based method; hence, we do not pursue its development further. 


E. SUMMARY 

This chapter concludes our development of NRBCs for the linearized Euler 
equations. For each NRBC, we have derived its implementation and demonstrated 
its effectiveness for short time-integrations. In the next two chapters, we will consider 
the question of longer time-integrations and the relative strengths and weaknesses of 
each method. 
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J 

E, 

Eu 


Ep 

1 

0.059866 

0.16273 

0.038642 

0.059866 

2 

0.059442 

0.16184 

0.038079 

0.059443 

3 

0.059306 

0.16156 

0.037884 

0.059306 

4 

0.059252 

0.16144 

0.037804 

0.059252 

5 

0.059236 

0.16139 

0.037775 

0.059237 

6 

0.059234 

0.16138 

0.037765 

0.059235 

7 

0.059235 

0.16138 

0.037761 

0.059236 

8 

0.059236 

0.16138 

0.03776 

0.059237 

9 

0.059237 

0.16138 

0.037759 

0.059237 

10 

0.059237 

0.16138 

0.037759 

0.059238 


Table XXXV. Error norms (IV. 14) for J G 1... 10, wave-like Hagstrom-Warburton 
NRBCs, in fini te channel, basic system with no advection 


J 

E, 

Eu 

Eu 

Ep 

1 

0.053965 

0.078277 

0.040606 

0.053965 

2 

0.053561 

0.077765 

0.040013 

0.053561 

3 

0.05347 

0.077625 

0.039873 

0.05347 

4 

0.053457 

0.077589 

0.039842 

0.053457 

5 

0.053456 

0.077575 

0.039836 

0.053456 

6 

0.053455 

0.077568 

0.039835 

0.053455 

7 

0.053454 

0.077564 

0.039835 

0.053454 

8 

0.053453 

0.077562 

0.039835 

0.053453 

9 

0.053453 

0.077561 

0.039834 

0.053453 

10 

0.053453 

0.077561 

0.039834 

0.053453 


Table XXXVI. Error norms (IV. 14) for J G 1... 10, wave-like Hagstrom-Warbnrton 
NRBCs, infinite channel, basic system with left-to-right advection 
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VII. 


LONG-TIME STABILITY 


A. OBSERVATIONS 

In this chapter we address the stability of the NRBCs for long time-integrations. 
This issue is critical to numerical weather prediction and other applications which re¬ 
quire time frames beyond that of a single wave’s propagation through the domain. 

Higdon [68] discusses the numerical stability of his NRBC formulation, us¬ 
ing criteria developed by Kreiss [81] and Gustafsson et ah [52] (and Higdon’s own 
characteristic-based interpretation thereof [63]). While his scheme meets the de- 
hned stability criteria, there is still concern that a scheme which is stable for the 
single-variable Klein-Gordon equation is also stable long-term for the equivalent 
linearized Euler system. Long-term stability, surprisingly, has not been much dis¬ 
cussed or demonstrated in NRBG development. In the Givoli-Neta papers explor¬ 
ing the automated Higdon NRBGs and the Givoli-Neta auxiliary variable NRBGs 
[39, 40, 41, 42, 43, 90, 113, 114, 115, 116], only [90] discusses the long-time stability 
of the solution, showing that the automated Higdon NRBG of order J = 10 is stable 
for long time-integrations. In the papers exploring the Hagstrom-Warburton auxil¬ 
iary variable scheme [57, 37, 53, 55], long-time stability was claimed in [57] (although 
the numerical example plotted showed an error norm that was increasing over time), 
and [55] uses an evanescent mode correction (a second set of auxiliary variables) to 
ensure long-time stability. 

As for the papers which discuss various NRBG techniques for the linearized 
Euler equations, only the absorbing layer methods [1, 17, 72, 73, 88] discuss or demon¬ 
strate long-time stability. Specihcally, the numerical examples in [17] are carried out 
for long time-integrations, and the PMLs [1, 72, 73] require additional terms to make 
them stable ([88] claims medium-term stability and implies stability over long time- 
integrations). 

All three types of NRBGs presented in this dissertation suffer from some form 
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of long-time instability. One example may be seen in Fig. 44, which is the state 
variable u from the nnmerical example of Sec. IV.B.3, with J = 10, rnn nntil t = 100 s. 
As we can see, the valnes on the bonndary have become catastrophically large, and 
they have pollnted the domain, completely overwhelming the actnal wave (faintly 
visible near the center of the domain). Tables XXXVII-XXXIX list the maximnm 
order J for which the given NRBC is stable for the varions domain conhgnrations 
for short (24 s), mediu m (100 s), and long (10,000 s) time-integrations. Table entries 
marked “No” indicate that the conhgnration was nnstable even for J = 1. The 
npper limit to J for the short-term stability mirrors a condition seen while prodncing 
the resnlts given in [90]; even thongh the Higdon scheme is theoretically stable in the 
discrete case, ronnd-off errors in the hnite-precision implementation lead to instability. 
This same problem afflicted the examples given in Chapter IV: If the example was 
performed on a 50 x 50 grid, it was stable only np to J = 5; halving the grid spacings 
enabled the stable resnlts np to J = 10. 


X'direction Vetocity pefturtatron 



Fignre 44. Plot of nnstable u in basic system (IV.1) with Higdon NRBC J = 10 in 
a semi-inhnite channel integrated np to t = 100 s. Notice the faint wave crests near 
the middle of the domain 
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Domain 

Coriolis 


Short 

Medium 

Long 

shape 

or gravity 

Advection 

(24 s) 

(100 s) 

(10,000 s) 

Bucket 

None 

None 

10 

8 

4 


Coriolis 

None 

10 

8 

4 


Gravity 

None 

10 

8 

No 

Channel 

None 

L-R 

10 

4 

No 


Coriolis 

L-R 

2 

2 

No 


Gravity 

None 

10 

8 

4 



L-R 

10 

4 

No 

Open 

None 

None 

10 

5 

2 



L-R 

10 

4 

No 



BL-TR 

10 

4 

No 


Coriolis 

None 

10 

5 

2 



L-R 

2 

No 

No 



BL-TR 

2 

No 

No 

Ground 

Gravity 

None 

10 

5 

No 



L-R 

10 

4 

No 


Table XXXVII. Higdon NRBCs, maximum stable order J for various domain config¬ 
urations and simulation durations 


For the Hagstrom-Warburton NRBCs, we found short-term stability as high 
as J = 40. While the NRBC still exhibits exponential improvement (see Fig. 45), 
the improvement is slight enough that the computational overhead for such high J is 
unjustihed; hence, we did not seek an actual “maximum” stable J but simply noted 
that the maximum is at least J = 40. 
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Domaiu 

Coriolis 


Short 

Medium 

Loug 

shape 

or gravity 

Advectiou 

(24 s) 

(100 s) 

(10,000 s) 

Bucket 

Noue 

Noue 

15 

7 

No 


Coriolis 

Noue 

15 

7 

No 


Gravity 

Noue 

13 

6 

No 

Chauuel 

Noue 

L-R 

12 

7 

4 


Coriolis 

L-R 

11 

7 

1 


Gravity 

Noue 

15 

7 

No 



L-R 

12 

7 

3 

Grouud 

Noue 

Noue 

15 

7 

No 



L-R 

12 

3 

No 


Coriolis 

Noue 

15 

7 

No 



L-R 

10 

3 

No 


Gravity 

Noue 

13 

No 

No 



L-R 

11 

3 

No 


Table XXXVIII. Givoli-Neta NRBCs, maximum stable order J for various domaiu 
coufiguratious aud simulatiou duratious 


Domaiu 

Coriolis 


Short 

Medium 

Loug 

shape 

or gravity 

Advectiou 

(24 s) 

(100 s) 

(10,000 s) 

Bucket 

Noue 

Noue 

40+ 

40+ 

No 


Coriolis 

Noue 

40+ 

40+ 

No 

Chauuel 

Noue 

L-R 

40+ 

40+ 

40+ 


Coriolis 

L-R 

40+ 

40+ 

40+ 


Gravity 

Noue 

40+ 

40+ 

No 



L-R 

40+ 

40+ 

5 

Grouud 

Noue 

Noue 

40+ 

28 

No 



L-R 

40+ 

3 

No 


Coriolis 

Noue 

40+ 

27 

No 



L-R 

40+ 

4 

No 


Table XXXIX. Hagstrom-Warburtou NRBCs, maximum stable order J for various 
domaiu coufiguratious aud simulatiou duratious 
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Density perturbation error norm 


X-velocity perturbation error norm 



Y-velocity perturbation error norm Pressure perturbation error norm 



Figure 45. Logarithmic plot of error norms (IV. 14) for J G 1... 40, Hagstrom- 
Warburton NRBCs, open half-plane with Coriolis, no advection. (TL) Error norms 
for p. (TR) Error norms for u. (BL) Error norms for v. (BR) Error norms for p. 


135 






B. SPECULATIONS 

Based on observations from these and other experiments dnring this research, 
we can tentatively identify several possible instability sonrces. These sonrces, thongh 
cnrrently mere conjectnre and specnlation, will be explored in fntnre research. 

1. Round-off errors. Early tests of these NRBCs nsed a smaller domain, only 1 
km sqnare, and the scheme was stable np to J ~ 8. Increasing the domain 
size and grid spacing by an order of magnitnde made the scheme stable np to 
J = 10. 

2. Normal derivatives reaching too far into the domain. The hrst nnmerical ex¬ 
amples, nsing the Higdon scheme, were rnn on a 50 x 50 grid. The scheme was 
nnstable for J > 5. Changing the grid to 100 x 100 made the scheme stable 
np to J = 10. As a resnlt, the normal derivatives of the Higdon scheme did 
not reach as far into the domain. Van Joolen [113] noted this same behavior 
for the Klein-Gordon eqnation. 

3. Dependency between interior and NRBC near corners. Althongh no points 
depend on the corner valnes in the Higdon NRBC scheme (see Fig. 13), the 
interior point closest to the corner depends on two bonndary points, rather 
than jnst one. It conld be that this dependence creates a slight error, jnst 
enongh that it grows over time and destabilizes the system. 

4. Frequency mis-match between interior and Higdon NRBC schemes. A reviewer 
for [21] noted that the Higdon NRBC discretization scheme (IV.10) cannot re¬ 
solve the shortest wavelengths resolvable by the interior scheme (IV. 7). When 
a wave strikes the open bonndary, the NRBC compntes the valne based on 
every other point at every other time step, as if for a wave with twice the 
wavelength. The skipped points of the wave are then nsed at the next time 
step to compnte that next time step’s bonndary valne. Hence, a wave striking 
the bonndary with wavelength n is resolved by the NRBC as a two-stage com¬ 
position of two waves, each with wavelength 2z/. A stability analysis by one of 
the co-anthors of [21] is in progress to determine the impact of this resolntion 
discrepancy. 

5. Numerically-singular matrix computations. For the Givoli-Neta NRBGs with 
gravity, in either the “bncket” or open-air domains, Matlab issned “Matrix is 
close to singnlar or badly scaled” warnings for large J and small St. It appears 
the left-hand matrix nsed to solve the anxiliary variables is prone to nnmeric 
instability. 

6. Exponentially-growing solutions. Althongh Higdon [68] proves the stability 
of his system by demonstrating that exponentially-growing solntions are not 
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permitted for the discrete system and boundary condition, it is possible that 
such solutions are appearing in this system. Look at Fig. 46, which shows the 
results of an unstable conhguration run until t = 100s. For each unstable J, the 
cx)-norm of the solution grows exponentially after a certain point in time. If we 
were to extend these lines backward, they would intersect at approximately 
t = 0. Hence, it appears that there is an initially small, but exponentially 
growing, solution within the scheme. Over time, this exponentially growing 
solution overwhelms the other interior values and leads to the unstable system. 

7. Discrete reflection coefficient greater than unity. Our analysis of the reflection 
coefficients (III. 15) and (III.34) is based on the partial derivatives in the con¬ 
tinuous case. However, we implement the system using hnite differences on 
discrete variables. It is possible that the “reflection coefficient” in the discrete 
case is dependent on more factors than just the wave speed and the NRBC 
estimated wave speed, factors such as the grid spacing and the time step size. 
Some conhgurations may lead to ||i?|| > 1 for some waves. See Appendix D 
for some initial analysis. A full analysis of this discrete reflection coefficient 
will be performed in future research. 



Figure 46. Logarithmic plot of state variable oo-norms for J G 1... 11, Givoli-Neta 
NRBCs, open half-plane, left-to-right advection, no gravity or Coriolis, run until 
t = 100s. (TL) cx)-norms for p. (TR) oo-norms for u. (BL) oo-norms for w. (BR) 
oo-norms for p. 


One result of the Hagstrom-Warburton table was surprising: The “bucket” 


conhguration was always unstable long-term, but the inhnite channel conhguration 
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showed long-term stability. This result was not a programming error. Running the 
simulation on a horizontal semi-inhnite channel showed the same instability as the 
vertical semi-inhnite channel. The stability of the boundary condition was impacted 
by the type of boundary on the opposite side. Until a formal analysis is undertaken 
to explain this result, we offer here a tentative hypothesis. The instability comes 
from waves whose rehection coefficients are slightly greater than unity when striking 
an open boundary on one side. On an open boundary on the other side, the wave 
speed and NRBC wave speed combine to cause a rehection coefficient less than unity; 
however, if the other side is a hard wall, then its rehection coefficient is exactly unity, 
and so the undiminished wave is totally rehected back to the open boundary, where 
it is again rehected and magnihed back toward the hard wall, and so forth, growing 
in magnitude with each pass. 


138 



VIII. 


SUMMARY AND COMPARISONS 


Having developed three distinct NRBC methods, we now examine them side- 
by-side to assess their relative strengths and weaknesses. As stated in the introdnction 
(Chapter I), there are fonr criteria we desire for an NRBC: speed, accnracy, stability, 
and ease of implementation. We consider each criterion individnally and compare the 
three NRBC methods against each other. 

Speed. We do not have speed comparisons of the implementations, as a simple 
“execntion time” metric is dependent on the efficiency of the code and the inventive¬ 
ness of the programmer. However, we can make some ballpark estimates based on 
the approximate operator connts reqnired. As noted in their respective chapters, the 
basic Higdon scheme (as antomated by the Givoli-Neta algorithm) reqnires 0(3'^) 
operations, bnt the simplihcation of setting all the cj to the same valne rednces the 
operation connt to O(J^). The Givoli-Neta and Hagstrom-Warbnrton anxiliary vari¬ 
able methods each reqnire 0{J) operations, except for the Givoli-Neta NRBC in the 
presence of gravity, which reqnires O(J^). 

Accuracy. By nsing the same nnmerical example thronghont this disser¬ 
tation, we can easily compare the relative accnracy of the three implementations. 
The only change is that the Higdon scheme’s no-advection examples nsed a semi- 
inhnite channel instead of the inf ini te channel nsed for the Givoli-Neta and Hagstrom- 
Warbnrton schemes. Even if we assnme the error norms for the Higdon scheme wonld 
be donbled by having two open bonndaries, we still see that the Higdon NRBC error 
norms are approximately 80% lower than Hagstrom-Warbnrton NRBC error norms, 
which in tnrn are approximately 40% lower than the Givoli-Neta NRBC error norms 
(for the basic system, bnt only 10% lower for the Coriolis- and gravity-inflnenced 
systems). The only exception to this case is the Higdon scheme with advection and 
Coriolis. In all other cases, Hagstrom-Warbnrton is approximately the same or slightly 
better than Givoli-Neta, bnt Higdon is signihcanly better than either of them. Again, 
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this is a consequence of the characteristic-based boundary method. 

Stability. When we consider stability, our comparison depends on the con- 
hguration. Comparing the Higdon scheme to the Givoli-Neta scheme, we see that 
Higdon is more stable when there is no advection, but Givoli-Neta is more stable in 
the presence of non-zero advection. In fact, the Higdon scheme is completely unstable 
in the non-zero advection case with Coriolis forces; however, the Givoli-Neta scheme 
shows some medium-term stability. The Hagstrom-Warburton scheme is the most 
stable, although it and the Givoli-Neta schemes are less stable in an open half-plane 
with non-zero advection. It also turns out that only the Hagstrom-Warburton scheme 
exhibited any appreciable long-term stability (specihcally, in the inhnite channel do¬ 
main), while the other two methods failed for nearly all conhgurations. 

Ease of implementation. All three methods are fairly straightforward and 
easy to implement. Since all three methods are explicit, they can be computed in 
a separate calculation after the interior scheme’s calculations. The Higdon scheme 
requires a summation algorithm to automate the high-order hnite-difference approx¬ 
imations, and the Givoli-Neta and Hagstrom-Warburton schemes require a careful 
consideration of characteristic and some matrix calculations. These drawbacks are 
not signihcant. However, Higdon is limited by the number of points in the domain, and 
it is not yet possible to implement Hagstrom-Warburton on an open top boundary in 
the presence of gravity. Furthermore, Higdon requires a careful choice of interior dis¬ 
cretization schemes [39, 19], while experiments with the Hagstrom-Warburton scheme 
for the scalar wave equation [55] show that higher-order interior schemes can be used 
(an option we have not attempted to utilize here). Such high-order schemes may also 
be permissible for the Givoli-Neta method. 

Summary. The Higdon scheme has an enormous advantage in accuracy, 
Hagstrom-Warburton is the most stable, Givoli-Neta and Hagstrom-Warburton are 
approximately equally fast, and all three are equally easy to implement. 
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IX. CONCLUSIONS AND AREAS FOR 
FURTHER RESEARCH 

We have shown through the preceding chapters that high-order non-reflecting 
boundary conditions can be applied to the 2-D linearized Euler equations in a wide 
variety of configurations. The Higdon, Givoli-Neta, and Hagstrom-Warburton tech¬ 
niques can all be used, in channels and in open domains, with and without advec- 
tion, with basic conditions or with gravity or Coriolis, excepting only the Hagstrom- 
Warburton method with gravity. However, despite the large amount of work con¬ 
tained herein, there remain many more extensions and improvements. These areas 
for further study include the following: 

1. Application of the auxiliary variable methods (Givoli-Neta and Hagstrom- 
Warburton) to finite element models 

2. Extension of the NRBC formulations to account for evanescent modes as well 
as the primary waves 

3. Analysis of the stability of each scheme for long-time integrations, including 
an exploration of the conjectures enumerated in Chapter VH, and mitigation 
of identified error sources 

4. Application of the NRBCs to the full three-dimensional system, including both 
gravity and Coriolis simultaneously 

5. Incorporating the NRBCs into a nested environment 

6. Finding a means to incorporate gravity into the Hagstrom-Warburton scheme 

7. Application of the NRBCs to the non-linear system, in two or three dimensions 

8. Extending this work to other linear first-order systems, such as Maxwell’s 
equations or the shallow-water equations (as a first-order system, not converted 
to the Klein-Gordon equation as in [113] and elsewhere) 

In the introduction, we stated that one motivation behind this NRBC develop¬ 
ment was to support the next generation of atmospheric modeling tools. This research 
has made significant progress toward that goal. With a broad-based finite-difference 
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implementation, the next step is to adapt the implementation for a spectral element 
system. In addition, the long-term stability concerns still need to be addressed. Al¬ 
though we have not satished our original purpose, we have nonetheless developed new 
computation tools for a broad array wave propagation models. 
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APPENDIX A. SIMPLIFYING THE EULER 

EQUATIONS 


1. INTRODUCTION 

Eq. (11.55) from Sec. II.B is 

dtp + dx{pu) + dy{pv) + d^{pw) = 0 

dt{pu) + d^{pu^) + dy{puv) + d:,{puw) + d^p = fpv 
dt{pv) + dxipuv) + dy{pv^) + d^{pvw) + dyp = -fpu (A.l) 
dt(pw) + dx{puw) + dy{pvw) + dz{pw^) + d^p = -gp 
dt{pe) + d^{{pe + p)u) + dy{{pe + p)v) + d^{{pe + p)w) = -gpw 


This form can be simplified, if we assume the fluid under consideration is an 
ideal gas {p = pRT). The simplihed form requires fewer terms and is thus easier and 
faster to calculate. We will consider the mass, momentum, and energy equations 
separately. 

2. MASS EQUATION 

We begin the simplihcation process with the hrst equation of the set: 

dtp + d^{pu) + dy{pv) + d^{pw) = 0 (A.2) 

This equation needs no simplihcation. However, we will expand the three spatial 
derivative terms using the product rule, since that form will appear in the subsequent 
simplihcations 

dtp + udxp + vdyp + wdzp + p {dxU + dyV + d^w) = 0 (A.3) 

3. MOMENTUM EQUATIONS 

The next simplihcation can be taken with the next three equations of the set: 
dt{pu) + d^{pu^) + dy{puv) + d^{puw) + d^p = fpv 
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} = fpv 

fpu 


> = 


> = 


-9P 


(A.5) 


dt{pv) + d^{puv) + dy{pv^) + d^{pvw) + dyp = -fpu 
dt{pw) + da;{puw) + dy{pvw) + d^{pw‘^) + d^p = -gp (A.4) 

We first use the product rule to separate the equations’ components as follows: 

udtp + udx{pu) + udy{pv) + udz{pw) 

+pdtu + pudxU + pvdyU + pwdzU + dxP 

vdtp + vdx{pu) + vdy{pv) + vdz{pw) 

+pdtv + pudxV + pvdyV + pwdzV + dyp 

wdtp + wdx{pu) + wdy{pv) + wdz{pw) 

+pdtw + pudxW + pvdyW + pwdzW + dzP 

We then combine terms to get 

u {dtp + dx{pu) + dy{pv) + dz{pw)) 

+p {dtu + udxU + vdyU + wdzu) + dxP 

V {dtp + dx{pu) + dy{pv) + dz{pw)) 

+p {dtv + udxV + vdyV + wdzv) + dyP 

w {dtp + dx{pu) + dy{pv) + dz{pw)) 

+p {dtw + udxW + vdyW + wdzw) + dzP 
The terms in the hrst set of parentheses matches (A.2), so we can eliminate them. 
Dividing the remaining terms by p gives 

dtU + udxU + vdyU + wdzU H —dxP = fv 

P 

dtV + udxV + vdyV + wdzV ^—dyp = —fu 

p 

dtW + udxW + vdyW + wdzW H —dzP = —g 

P 


I = fpv 
\ = -fpu 
= -gp 


(A.6) 


(A.7) 


4. ENERGY EQUATION 

Now the simplihcation process gets more complicated, using the hnal equation 
of the set: 

dt{pe) + dx{{pe + p)u) + dy{{pe + p)v) + dz{{pe + p)w) = -gpw (A.8) 
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We begin with our state equation and our definition of energy 


p = pRT 

(A.9) 

e = c^T + ^{u^ + + w"^) 

Li 

(A. 10) 

R = Cp — Cy 

(A.11) 

Combining these terms, we have 


Cv pu^ + pn^ + ptc^ 

2 

(A. 12) 

Substituting this value into (A.8) gives 



dt (pf + 

+d, ((pt + u + pu) 

+dy ((pf + £!i!±£|!W) V + pv^ 

+d, ((pf + e^!±£|!W) ^ + pui^ 

Separating the sums and using the product rule gives the expansion 

^dtp + pudfU + pvdtv + pwdtw + 

+ %dx{pu) + (pw + W^)) + 5a:(pw) 

+^ 5 j^(pn) + \dy {pv {u^ + + w^)) + dy{pv) 

+%dz{pw) + {pw {v? + + w^)) + dzipw) 

Combining the ^ terms and expanding more pieces with the product rule gives 

\ 

^ {dtp + udxP + vdyp + wdzp + p {dxU + dyV + d^w)) 

+pudtu + pvdtV + pwdtW 
+ ^2+,P^,^2 aa;(pM) + 5j^(pn) + 5 ^(pm;)) 

+(pm) (u^a^M + vdxV + wdxw) ' = —Qpw (A. 15) 

+ {pv) {udyU + VdyV + WdyW) 

+{pw) {udzU + + wd^w) 

+udxP + vdyp + m;5^p + p (5a;M + 5^^ + d^w) 
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The parenthetical term in the third line matches (A.2), so we can eliminate it entirely. 
In addition, if we combine the pu, pv, and pw terms together, we get 

^ {dtp + udxP + vdyp + wdzP + p {dxU + dyV + dzw)) 

+{pu) {dtu + udxU + vd^v + wdxw) 

+{pv) {dtv + udyU + vdyV + wdyw) * = —gpw (A. 16) 
+{pw) {dtw + udzU + vdzV + wdzw) 

+udxP + vdyP + wdzP + p {dxU + dyV + dzw) 

Based on the simplihed resnlts in (A. 7), we can replace the third, fonrth, and hfth 
lines of the above eqnation, resnlting in the following 

^ {dtp + udxP + vdyp + wdzP + P {dxU + dyV + dzw)) 

+W {-),dxP + fv) 

+(pn) (-),dyp - fu) > = -gpw, (A. 17) 
Hpw) (-),dzp-g) 

+udxP + vdyp + wdzP + P {dxU + dyV + dzw) 

which we can expand into 

^ {dtp + udxP + vdyp + wdzP + P {dxU + dyV + dzw)) 

—udxP + fpuv — vdyP — fpuv — wdzP — gpw = —gpw, (A. 18) 
+udxP + vdyP + wdzP + P {dxU + dyV + dzw) 

Canceling terms rednces this eqnation to 

{dtp + udxP + vdyP + wdzP + P {dxU + dyV + dzw)) + p {dxU + dyV + dzw) = 0, 

R 

(A.19) 

which we can also write as 


^ {dtp + udxP + vdyp + wdzp) + (^ + l) T + dyV + dzw) 


From onr dehnition of R we have 


C|; R Cy T Cy 

R^ ^ R ^ R 


(A.20) 

(A.21) 
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so that multiplying the above equation by — gives the simplihed equation 

C-v 

dtp + ud:cP + vdyp + wdzP + 7P {d^u + dyV + = 0, (A.22) 

where 7 = ^, and we have replaced our energy variable with the primitive state 
variable for pressure. 
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APPENDIX B. THE FINITE DIFFERENCE 
INTERIOR SCHEME 


Let us briefly analyze the interior scheme used in Chapters IV~VI. We claim 
that this discretization scheme is 0{Sx‘^ + Is that claim true? To test this 

scheme, we contrive an analytic solution of sines and cosines 


Po cos{kxx) cos{kyy) cos(u;t) 
Mo sin(/ca;a;) cos{kyy) cos(a;t) 
Vo cos(kxx) sin(kyy) cos(a;t) 
Po cos(kxx) cos(kyy) cos(ujt) 


(B.l) 


with kx = ky = j and oo = • \Jkx + ky. On a 40 x 40 m domain, these values 

match a hard wall boundary condition on all four sides. We apply the differential 
operators of (IV. 1) to this equation, and set the results as the right-hand-side values of 
(IV.1). These terms will then act to force the solution to remain equal to (B.l) within 
the limits of the discretization error. So that all four state variables are approximately 
the same order of magnitude, we set po = 1; 7 = 1; Po = 1; and Mq = = 0.25. 

Since the leap-frog method requires two prior time steps, we use the analytic 
solution to set the values for the hrst two time steps. We then use the discretization 
scheme to compute the next time step, which we compare to the analytic solution at 
that same time step. We compute the average absolute error at each interior point 


for each state variable by 


(N, - 2 ) (N, - 2 ) 


E^p — 


(B.2) 


where p is our computed state variable, and (p is the analytic solution state variable. 

For the hrst test, we begin with 5x = 0.8, and we halve it with each iteration. 
We set 5y = 0.00625 for all iterations, so that the discretization error in y does not 
overwhelm the ^-discretization error we wish to measure, and we set 5t = 2“®, which 
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Sx 

Sy 

St 

Ep, Ep 

Eu 

Ex 

0.8 

0.00625 

0.015625 

2.7983 xl0~^ 

8.0102 xl0~^ 

4.9978 xlO-^ 

0.4 

0.00625 

0.015625 

7.1 xlO"^ 

2.0223 xlO-^ 

1.2617 xlO"® 

0.2 

0.00625 

0.015625 

1.7814 xlO-5 

5.062 xlO-^ 

3.1584 xlO-® 

0.1 

0.00625 

0.015625 

4.4562 xlO-6 

1.2654 xlO-® 

7.91 xlO'^ 


Table XL. Discretization errors for different grid spacings Sx 


Sx 

Sy 

St 

Ep, Ep 

Eu 

Eu 

0.00625 

0.8 

0.015625 

2.7983 xlO-^ 

4.9978 xl0~^ 

8.0102 xlO-4 

0.00625 

0.4 

0.015625 

7.1 xlO-^ 

1.2617 xlO-® 

2.0223 xlO-^ 

0.00625 

0.2 

0.015625 

1.7814 xlO-® 

3.1584 xlO-6 

5.062 xlO-5 

0.00625 

0.1 

0.015625 

4.4562 xlO-6 

7.91 xlO-^ 

1.2654 xlO-® 


Table XLI. Discretization errors for different grid spacings Sy 

is well below the CFL limit for these conditions. Table XL shows the discretization 
errors for each state variable with each halving of Sx. The error norms decrease by 
a factor of almost fonr with each halving of Sx, exacly as expected for a second- 
order method. The only exception to this decrease is the error norm for v. We note, 
however, that the eqnation for v contains two dy terms and only one dx term, so it 
is not snrprising that redncing the error for dx without simultaneously reducing the 
error for dy wonld have only a small impact. 

For onr second test, we set Sx = 0.00625, and we begin with Sy = 0.8, halving 
it with each iteration. Keeping St the same as before, we get the resnlts shown in 
Table XLI. Note that these errors are almost identical to those for the Sx test, with 
the error norms for u and v reversed. Hence we conclnde that the method is also 
second-order in the y direction. 

As a third test, we start with Sx = Sy = 0.8 and halve both grid spacings at 
each iteration. Using a larger St = 2~^, we get the resnlts shown in Table XLII, and 
we see second-order resnlts in all fonr state variables. 

Finally, we test St. This time, we set Sx = Sy = 0.4, and we begin with 
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Sx 

Sy 

St 

Ep, Ep 

Eu 

E^ 

0.8 

0.8 

0.03125 

9.5561 xl0~^ 

1.5883 xlO-^ 

1.5883 xlO-^ 

0.4 

0.4 

0.03125 

2.4401 xlO-^ 

4.031 xlO-^ 

4.031 xlO-^ 

0.2 

0.2 

0.03125 

6.1292 xlO-^ 

1.012 xlO-^ 

1.012 xlO-^ 

0.1 

0.1 

0.03125 

1.5216 xlO-5 

2.5368 xlO-® 

2.5368 xlO-® 


Table XLII. Discretization errors for different grid spacings Sx and Sy 


Sx 

Sy 

St 

Ep, Ep 

Eu 

Eu 

0.8 

0.8 

0.25 

0.0066075 

0.012469 

0.012469 

0.8 

0.8 

0.125 

0.0037399 

0.0063109 

0.0063109 

0.8 

0.8 

0.0625 

0.0019051 

0.0031717 

0.0031717 

0.8 

0.8 

0.03125 

0.00095561 

0.0015883 

0.0015883 


Table XLIII. Discretization errors for different time steps St 


St = 0.25, halving it with each iteration. The resnlts are somewhat snrprising. As 
shown in Table XLIII, the error norms only decrease by a factor of two with each 
halving of St, which implies a hrst-order method rather than second-order. The 
discretization scheme is the same. How can a scheme which is second-order in space 
be only hrst-order in time? Perhaps the difference comes from how the discretization 
scheme is nsed. For the spatial derivative approximations, we nse the known node 
valnes to approximate the derivative, that is. 


dxU ^ 


u. 


■i+l 


u: 


-1 


2Sx 


(B.3) 


which is a second-order approximation, easily demonstrated via Taylor series expan¬ 
sion. For the time derivative, we approximate it nsing the eqnation system, and then 
nse that approximation and the earlier node valne to compnte the new node valne: 




2St 


dtu = RHS((p) 

i dtu 


,n+l 


+2(Si9,M 


(B.4) 

(B.5) 

(B.6) 
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If we rewrite the superscript for the ^ term, we see that this method is in fact 
Euler’s method over 25t: 

~ + 28tdtu (B.7) 

Since Euler’s method is only a hrst-order approximation, our time marching method 
is only hrst-order, even though it is dehned using the same discretization scheme 
for our second-order spatial derivative approximations. Thus, it appears that our 
discretization scheme is in fact 0{5x^ + Sif St). 

This analysis “feels” wrong, and it caused some strenuous debates with the 
dissertation supervisors. However, it hts the observations. Furthermore, other ex¬ 
periments designed to test the leap-frog scheme’s performance also showed similarly 
inexplicable results. Testing it as an ODE solver resulted in 0{Sx^) performance for 
a single-step and a global 0(5a:^) convergence. Contrariwise, a test using an analytic 
solution to the heat equation 


OtU 

initially showed 0{5x^f‘^) single-step convergence, but as the time step became smaller, 
the improvement shrank to 0{5x) and remained there. The cause of this performance 
is under investigation. 
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APPENDIX C. WAVE-LIKE SOLUTIONS OF 
THE LINEARIZED EULER EQUATIONS 


We now derive the wave-like solution used in Lemma VL2. The derivation is 
similar to Hu’s [71], where he uses a non-dimensionalized system and split variables 
to derive the PML equations. We hrst assume the existence of a wave-like solution, 
then we derive the amplitude of each state variable. By avoiding contradictions, we 
prove that such a wave-like solution exists. This wave-like solution has the form 

^^^^ikx+ily-iu^t ^ (C.l) 


where (p denotes our state variables, and denotes the amplitudes of each compo¬ 
nent. If we apply (C.l) to (IV.1) and cancel out the common exponential term, we 
get 

—iojp* + ikuop* + ilvop* + ikpou* + ilpov* = 0 
—iu)u* + ikuou* + UvqU* + ^p* = 0 

(C.2) 

—iu!v* + ikuov* + Uvqv* -|- ^p* = 0 
—iujp* -\- ikuop* + ilvop* + ik'jpou* + iljp^v* = 0 . 

Combining terms, we get 


{kuo + Ivq — ou) p* + kpou* + Ipov* = 0 
{kuo -f Ivq — u)u* + Vp* = 0 
{kuo + Ivo -uj)v* + Vp* = 0 
k'ypou* + l^pov* -I- {kuQ + Ivq — u})p* = 0 . 


For acoustic waves, kuo + Ivq — co 0. {kuQ + Ivq — u) = Q for entropy and vorticity 
waves; see [71]). We can easily solve for p* and p* in terms of u* and v*: 


Po {kul + Iv*) 
oj — kuo — Ivq 
7Po {kul + Iv*) 
OJ — kuQ — Ivq 


(C.4) 
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From (C.3b,c) we have 


kuo + lvo — CO , kuo + lvo — CO , 
-:- u = -:- V 


(C.5) 


which leads directly to 

u* _ k 

V* I 

If we insert our solution for p* into (C.3b,c), we get 

* kp* 


Po 

— kuo 

- 

■ ivo) 

k^Po{ku* 

+ 

Iv*) 

Po 

— kuo 

- 

■ Ivof 

IPo 

k‘^u* 

+ 

klv* 

Po (a; — kuo 

- Ivof 


Ip* 



Po 

— kuo 

- 

■ Ivo) 

Ijpo {ku* 

+ 

Iv*) 

Po 

— kuo 

- 

■ Ivof 

7Po 

klu* 

+ 

t^v* 


Po {co — kuo — Ivof 


(C.6) 


(C.7) 


(C.8) 


If we use (C.6) to remove I from the numerator of (C.7a) and k from the numerator 
of (C.7b), we get 


7Po + 

Po (co — kuQ — IvoY 

IPO 


u = 


V = 


Po (co — kuQ — Ivof ’ 
which, after a little algebra, can be reformulated as 


(C.9) 


_^co — kuo — Ivo u* 

Co + v*'^ 

_^co — kuo — Ivq V* 

Co y/u*'^ + 



(C.IO) 
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Now, let u* = U cos 9, v* = U sin 6, for some U & C and 9 G [0,27r] measnred 
connterclockwise from the positive x axis. After some manipnlation, we can explicitly 
solve for k and 1: 


U! COS 9 

Co + Uq cos 9 + Vq sin 9 
u) sin 9 

Co + Uq cos 6* + no sin 6* ’ 


(C.ll) 


where we choose the positive for each ± in (C.IO), denoting propagation in the positive 
X and y directions. We then insert these valnes of k and I into the eqnations for 
and simplify, yielding 

p* = 

u* = U cos 9 

(C. 

V* = U sin 9 

p* = 

^ Co 

Thns onr wave solntion becomes 



0 


if exp 



cos X + sin y 


Co + Uq cos 9 + Vq sin 9 



(C.13) 


with dehned by (C.12). A little tedions algebra shows that this solntion satishes 
(IV.1). If we insert this solntion into (VI.16), we get the wave speed 


= (co + Uq cos 9 + Vq sin 9)'^ . 


(C.14) 


Similarly, when we consider the opposite-sign terms of (C.IO), we again get a wave-like 
solntion, this time with 

cIj = {uq cos 9 + Vq sin 9 — cq)'^ . (C.15) 


Since we do not know in advance the propagation angle of these plane waves, for 
the pnrpose of dehning the NRBC wave speed, we assnme the angle is normal to 
the bonndary. With this assnmption, and the reqnirement that c^i > 0, we have the 
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following: 




— Co + Mo 

Tn 

Cyj 

= Co + Mo 

Tw 

Cyj 

o 

1 

o 

Ts 

Cyj 

o 

1 

o 


(C.16) 
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APPENDIX D. DISCRETE REFLECTION 
COEFFICIENT^A PRELIMINARY ANALYSIS 


Let us briefly look at this idea of quantifying the reflection coefficient for the 
discrete case (see Chapter VII). This analysis differs from Sec. 8.1.5 of [24] in that 
we consider the physical waves, not merely the high-frequency non-physical compu¬ 
tational waves generated by the hnite-difference scheme. For our initial analysis, we 
consider a Sommerfeld (hrst-order Higdon) boundary condition in a 1-D domain. A 
full analysis is outside the scope of this dissertation. We begin the work here to show 
its potential for future research. 


1. DERIVATION 

As with the analysis of the reflection coefficient for the continuous equation 
(Sec. III.B.l.a), we begin by considering a wave of the form 

n(a;,t) (D.l) 

traveling left to right at unknown speed Cx- (For a two-dimensional wave, we consider 
only the portion of the wave traveling parallel to the x axis.) If we use a Sommerfeld 
condition {dt + Codx) u = 0, then the computed solution for u will include a reflected 
wave of the correct magnitude to satisfy the boundary condition. Thus, 

u{x, t) = (D.2) 


Discretizing the Sommerfeld condition with a backward difference in x and t, we have 


u 


N 


U 


n—1 

N 


U 


+ Cq- 


N 


U 


N-l 


h 


= 0 , 


(D.3) 


where the subscripts N and N — 1 denote the point on the right boundary and the 
point immediately to its left, respectively; the superscripts n and n — 1 denote the 
current and previous time steps, respectively; k denotes the time step size; and h 
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denotes the spatial step size. Thus, at the right boundary at time step n, we have 
X = Nh and t = nk. Applying this discretization to our wave (D.2) at this point, we 
have 

^i{Nh-Cxnk) _j_ h+c^nk) _ ^i(Nh—Cx(n-l)k) _ J^^i{Nh+Cx{n-l)k) 

+ A _j_ j^^i{Nh+c^nk) _ ^i({N-l)h-Cxnk) _ J^^i{{N-l)h-c^nk)\ 


= 0, (D.4) 


where A = ^. After we cancel and combine terms, we solve for R and get 

1 - + A (l - 

Q2icxnk _ Q-iCxk (^•^) 

We are interested in the magnitude of R more than its actual value. Solving for this 
magnitude, we have 


R 


R 


1 - + A (l - 

1 — + A (1 — e“*^) 

^_ sm{cxk) _ 

2i (1 - + A (1 - e-*^)) 


1 - 


sm c. 


A) 


— icxk 

e 2 sm 


Cxk 

2 


T Ae 


-ih 

2 


sm 


(I) 


(D.6) 


2. IMPLICATIONS AND SPECULATION 

Looking at (D.6), we see that certain combinations of c^k and h could make 
the fractional term negative, resulting in a reflection coefficient R with magnitude 
greater than one. In fact, if we use the following constants, 

m 

Co = 343— 
s 

h = 100 m 
k = 0.18 s , 

which are close to those used in our numerical examples in this dissertation, we get 
several large reflection coefficients as Cx varies from 1 up to cq. The top half of 
Fig. 47 shows these reflection coefficients. Each value of R larger than one represents 
a potentially unstable wave mode. 
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Constants for mesoscale Euler equations 
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Figure 47. Discrete reflection coefficients for varying wave speeds Cx- The x-axis is 
the wave speed Cx'-, the y-axis is the magnitude of the reflection coefficient R computed 
by (D.6). (Top) Discrete reflection coefficients using constants approximately equal 
to those in this dissertation’s numerical examples for a mesoscale model. (Bottom) 
Discrete reflection coefficients using the same constants as the numerical example of 
[90] for the Klein-Gordon equation in a small-scale model. 

On the other hand, if we use the values from the numerical example in [90], 

^m 

s 

= 0.25 m 
= 0.125 s , 

then our discrete reflection coefficients are all less than one, even if Cx is twice as large 
as cq. These coefficients are plotted in the bottom half of Fig. 47. These plots imply 
the stability of the example in [90] and the instability of our examples. 

The problem appears to be the scale of the domain. Perhaps the introduction 
of a scaling factor can improve the stability. Future research will expand this analysis, 
to determine the general formula for the Higdon NRBC of order J, and also to con¬ 
sider the Givoli-Neta and Hagstrom-Warburton NRBGs, after converting the normal 
derivatives to tangential. Ideally, this analysis will uncover the choice (or combina- 


co 

h 

k 
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tion of choices) of Cj to choose to keep ||i?|| < 1 for all possible c^. In the meantime, 
this qnick test has nncovered a possible reason for the instabilities exhibited by these 
NRBCs. 
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